{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-23T14:25:36.356812Z",
     "start_time": "2021-09-23T14:25:36.187362Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[9.99999000e-01 1.00000000e-06]\n",
      " [9.99995792e-01 4.20754569e-06]\n",
      " [9.99982488e-01 1.75123961e-05]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA23klEQVR4nO3de5zUdd3//8eT5SziLmfltGDgAYFVEUUz8YxadOUXU6PLU0WoXWldlfojs8viysuszFNEilZunsvQLMkTmpgKiiYqKoSwKnLaRc4Hef3+eH8Gh2V25jO7+9nZ3Xndb7e5zXw+8zm8ZnZ2XvN+vz/v91tmhnPOueLVptABOOecKyxPBM45V+Q8ETjnXJHzROCcc0XOE4FzzhU5TwTOOVfkPBG4nCT9VdK5dTxXLskktU3gvNMkXdnYx007/npJg7M8v0TSCUmdPxtJT0n6ah3PDYhiL8nzmCbpU40TYePJ9lqb87FbE08EzYikT0uaI2mtpDWSnpV0WPTceZL+kbbtEkmboi+E5ZLukNQlibjM7BQz+20Sx85x3slm9qMEj9/FzBYDRO/fj5M6V7raf8t8mdnSKPaPs5yj6L4AJf1Q0p2FjqMl8kTQTEjqCjwM3Ah0A/oC/wNsybLb58ysC1ABHAxckXCYrplT0Cz+r5MoJbpkNIsPjANgKICZ3WVmH5vZJjObZWav5trRzJYDjxISQkaS9pJ0m6QPJL0n6cepqoXoF+qzkm6MSiNvSjo+bd+dvy4llUi6TtIqSYuB0+pxnl9IqpG0WNKR0fplklakV0HV/pUu6fOS5kv6SNIiSeMyvM7zJT2UtvyOpHvTlpdJqogem6RPSZoETAS+F5WwHko7ZIWkV6P35R5JHdOO9bXo+GskzZS0T7R+t+qy1Hso6QBgGjAmOldNXX8zYF9JL0Tn/rOkbpmOHx17qqRngY3A74GjgZuic9yUdswTJL0tqVrSzZIUHeNTkmZH51ol6Z66gpI0XtKC6G/4VPSaUs8tkXSZpFeBDZmSgaQTo8/Y2ig21Xr+AklvRDE+Kmlg2nO/jP6GH0maJ+noaP044P8Dzoxe8ytphxwYfe7WSZolqUe0T0dJd0paHb2WFyX1zvL3aL3MzG/N4AZ0BVYDvwVOAcpqPX8e8I+05SXACdHjfsC/gF9mOf6DwK+BPYBewAvA19OOvR34FtAOOBNYC3SLnn8K+Gr0eDLwJtCfUHJ5EjCgbR7nOR8oAX4MLAVuBjoAJwHrgC7R9ncAP44ej45iOpHwA6YvsH+G1zkYqIm22Rt4F3gv7blqoE20bMCnap+r1nv8ArBP9FrfACZHzx0HrAIOiWK/EXg6eq48/T3J8B7u8res4+/1FPAecFD0Xj4A3Jnp+NG2S4FhQNvob7jzfGnHNEKpsxQYAKwExkXP3QVMid63jsCn64hrKLAh+ju0A74HvAO0T3vP5hM+H50y7N8D+AiYEO3/regzkXpv/iM63gHRa/k+MCdt/y8D3aPn/htYDnSMnvth6j2q9T4uiuLuFC1fEz33deAhoDPh83go0LXQ3wWFuHmJoJkws4+ATxP+WX8DrIx+ZWb7hfKgpHXAMmAFcFWmjaJjnAJcamYbzGwF8AvgrLTNVgDXm9k2M7sHWEitX/uRL0bbLTOzNcBP8jzPv83sdgv12/cQvjCuNrMtZjYL2ApkatD8CjDDzP5uZjvM7D0ze7P2Rhbq/NcRSkfHEEpK70naP1p+xsx2ZHqf6nCDmb0fvdaH+KTUNTGK5yUz20KolhsjqTyPY+fyezN7zcw2AFcCX1TdDcR3mNkCM9tuZtuyHPMaM6sxs6WEJF4Rrd8GDAT2MbPNZlZXG8aZwF+iv8M24DrCF+yRadvcEH0+NmXY/1TgdTO7P9r/esKXecrXgZ+Y2Rtmth34X0KpbCCAmd1pZquj1/kzQhLeL8vrBbjdzN6K4rm31mvuTvgx8LGZzYv+D4uOJ4JmJPrwn2dm/Qi/BPch/KPU5T/MbE9gLLA/4ddWJgMJv74+iIrANYRf7b3StnnPop9JkXej89e2DyHxpG+Xz3k+THu8CcDMaq/L1Ojdn/DLLo7ZhPfkM9HjpwhJ4JhoOR/pX1Ib02Lbh7TXbmbrCSW6vnkeP5va73M76v4bL6tjfW11vZ7vEapoXoiqfS6oY//ar3tHdO70150tll0+P9FnLn37gcAv0z4/a6K4+gJI+u+o2mht9Pxe1P2epNT1mn9P+KFwt6T3JV0rqV2OY7VKngiaqejX7h2EhJBr29nRttfVsckyQqNzDzMrjW5dzWxY2jZ9U/XFkQHA+xmO9QHhSzl9u3zOU1/LgH1jbptKBEdHj2eTOxHkOwzv+4QvLQAk7UH4dfkeoeoEQpVDSp96nKv2+7yNUB2VSe1j5vV6zGy5mX3NzPYh/Cq/RZkvNa39uhXF+V7Mc+/y+UnbP2UZoSqxNO3WyczmRO0BlxFKpWVmVkqoLkx9bvN9zdvM7H/M7EBCieazwDn5HKO18ETQTEjaP/q10y9a7g+cDfwz5iGuB05MNYSmM7MPgFnAzyR1ldRG0r6SjknbrBfwTUntJJ1BqKN9JMN57o226yepDLg8z/PU123A+ZKOj47bN6ruyWQ2cCyhjroKeAYYR/iifrmOfT4ktCHE9YcongpJHQhVGM+b2RIzW0n4YvyyQuP6BeyaxD4E+klqn+McX5Z0oKTOwNXA/ZblktGGvB5JZ6Q+e4R2FAMynete4LTo79COUE+/BZgT81R/AYZJOj1qSP4muybJacAVkoZFce0VfR4B9iS0J6wE2kr6AaFtLeVDoFwxr5qSdKyk4VF120eERBv3/W1VPBE0H+uAw4HnJW0gJIDXCP9oOUVfPr8j1CVncg7QHnid8I9+P6ExNeV5YAjhF+dUYIKZrc5wnN8QitOvAC8Bf8zzPPViZi8QGpl/QfgVOJu0X6a1tn0LWE9IAKn2l8XAs1m+SG8DDoyqJB6MEc/jhPf6AcKv3H3ZtS3ka8B3CdVFw9j1i/IJYAGwXFJdv/AhVF3cQdQgSvjSjOuXwIToypsbYmx/GOGztx6YCVxiZv+uvZGZLSQ02N5I+Kx8jnAZ89Y4QZnZKuAM4BrCezMEeDbt+T8B/0eorvmI8D9wSvT0o8BfgbcI1VOb2bVa6b7ofrWkl2KE04fw+fyIcCHAbKAo+yFo12phV4wknUe4auPThY7FOdf0vETgnHNFzhOBc84VOa8acs65IuclAuecK3ItblCoHj16WHl5eaHDcM65FmXevHmrzKxnpudaXCIoLy9n7ty5hQ7DOedaFEnv1vWcVw0551yR80TgnHNFzhOBc84VOU8EzjlX5DwROOdckUssEUiaoTD14Gt1PC9JNyhM9feqpEOSiKNPH5Ca7tanT+6YCqayEsrLoU2bcH/RRdmXKytb9j7NKRbfp3nF0tL3qazc/f+7IRp7yrPUjTApyCHAa3U8fyphJEEBRxCG8M153EMPPdTyAYW/9e6dV8iN4847zQYONJPC/YUXmnXunF/g7dqZtW/fcvdpTrH4Ps0rlpa+T+fO4X88D8Bcs8zfqzmHmJC0L1BlZlskjQVGAL8zs5pcSSaatu9hM9ttchVJvwaeMrO7ouWFwFgLY9rXadSoUZZPP4JdplopoBxvc8NVVsKUKbB0KXTrBuvWwdZYIwM751qigQNhyZLYm0uaZ2ajMj0Xp2roAeDjaLai24BBhEk5Gqovu44lXkUd0/xJmiRprqS5K1eubIRTtzKVlTBpErz7bsg4q1d7EnCutVu6tNEOFScR7LAwifQXCJOWf4tGmGiET6aXS5fxd7OZTTezUWY2qmfPjD2km71E2xGmTIGNGxv5oM65Zm3AgNzbxBQnEWyTdDZwLvBwtK4xJniuYte5SvuReY7cVufDD3Nvk5e4vwzyrSdr1w7a55pNsRnv05xi8X2aVywtfZ/OnWHq1PyOk0WcRHA+MAaYamb/ljSIxpnObSZwTnT10BHA2lztA/XRu3djH7EZilNK6twZJk8O9YpSuL/wwuzLt98OM2a03H2aUyy+T/OKpaXvM306TJzYaF8hic1HIOkuYCzQgzCp9FVEJQkzmyZJwE2EScU3AuebWc5W4Hwbi5tSnz7xf+03+G1Pbxw2Cx+Q9IO2awddu8KaNaEIOXVqo35wnHMtS7bG4pyjj0r6Nxnq7s1scLb9zOzsHM8bcHGu87cky5fvupzYFUupxuH0doG2bf2L3zlXL3GGoU7PIB2BM4BuyYTjYsnUOLxtG3TpAqtWFSYm51yLlbONwMxWp93eM7PrgeOSD63lq6t9osHtFnU1Djfi5WTOueKRMxFIOiTtNkrSZGDPJoitxVu+/JOugFddFaqKqqp2r0LKW12XjTXi5WTOueIR56qhn6XdfgIcCnwxyaBao4kTQ0K4665GONjUqVBSsuu6Rr6czDlXPBK7aigpzfmqoVwOPzx0+H355QYeyAzKykK7wKZN3jjsnMupQUNMSNpL0s9TQzxI+pmkvRo/zNZv4kSYPx8WLGjggZYuhbVr4dprYceOMN6IJwHnXD3FqRqaAawjVAd9EfgIuD3JoFqrM88MNToNHkF2zpxwf+SRDY7JOefiJIJ9zewqM1sc3f4HyNqHwGXWuzeceCL84Q/hh3y9Pfcc7LEHDB/eaLE554pXnH4EmyR92sz+ASDpKGBTsmG1Tuk9j9Pbenv3zvNKojlzYPTo0InMOecaKM43yYXAb6N2AQFrgPOSDKq1qmv4ibwGoduwITQ0XH55Y4TknHO5E4GZzQdGSuoaLX+UdFAui7lz4eOPYcyYQkfinGsl6kwEkr5sZndK+nat9QCY2c8Tjs1lkmooPuKIwsbhnGs1spUI9ojuM/UiblmdD1qTOXNg//2he/dCR+KcayXqTARm9uvo4WNm9mz6c1GDsWtqZuGKofHjCx2Jc64ViXP56I0x17kcGjwI3dtvh/mIvf+Ac64RZWsjGAMcCfSs1U7QFSjJvJfLJnWJ6McfQ4cOcNlleQ4P9Nxz4d4TgXOuEWVrI2gPdIm2SW8n+AiYkGRQrV1JCfTrV49Ro+fMgdLS0EbgnHONJFsbwWxgtqQ7zOzdJoypKPTvD8uW5bnTnDnhaqE2cWr0nHMunjjfKLdKKk0tSCqT9GhyIRWHAQPyKBFUVoYdXnstVA81eLAi55z7RJxE0MPMalILZlYN9EosoiLRv3+YpCbnmEOp+YlTxYe1a8OyJwPnXCOJkwh2SNo59ZWkgXg/ggYbMCBMJ5BzeIlM8xNv3BjWO+dcI4gz1tAU4B+SZkfLnwEmJRdScUjNKrlsGey9d5YNfX5i51zC4kxe/zfgEOAe4F7gUDPzNoIG6t8/3Of8Pvf5iZ1zCYszQ5mAccAhZvYQ0FnS6MQja+VS3+M5E8HUqWE+4nQ+P7FzrhHFaSO4BRgDnB0trwNuTiyiIlFaGuaWyXkJ6cSJMH36J3MPDBwYln1qSudcI4mTCA43s4uBzbDzqqH2iUZVBKQ8LiGdOBG6doWLL/b5iZ1zjS5OItgmqYToSiFJPYGGTLToIgMGxOxUtm0brFkDvfyqXedc44uTCG4A/gT0kjQV+Afwv4lGVST6949ZIli1Ktx7InDOJSDboHODzOzfZlYpaR5wPGGqyv8wszeaLMJWbMCA0I9gy5YwCF2dVqwI954InHMJyNaP4H7gUEmPm9nxwJtNFFPRSF1CWlUF++6bZUNPBM65BGVLBG0kXQUMrT1dJfhUlY0h/RJSTwTOuULJ1kZwFuFKodQw1LVvOUkaJ2mhpHckXZ7h+b0kPSTpFUkLJJ2f/0touVIlgpwNxp4InHMJyjYM9ULg/yS9amZ/zffA0ZVGNwMnAlXAi5JmmtnraZtdDLxuZp+LrkZaKKnSzLbme76WqF+/cJ+zwXjFCmjXDvbaK/GYnHPFJ85YQ09I+hJQnr69mV2dY7/RwDtmthhA0t3A54H0RGDAnlHv5S7AGmB77OhbuE6dwo/8nCWClSuhZ8/Q+cA55xpZnETwZ2AtMA/Yksex+wLpX3FVwOG1trkJmAm8T6huOtPMduujIGkS0UB3A1rZGDuxLiFdscKrhZxziYmTCPqZ2bh6HDvTz9faw1efDMwHjgP2Bf4u6Rkz+2iXncymA9MBRo0a1aqGwB4wAN56K8dGngiccwmK06FsjqTh9Th2FdA/bbkf4Zd/uvOBP1rwDvBvoKgm5O3fH959FyxbevNE4JxLUJxE8GlgXnT1z6uS/iXp1Rj7vQgMkTRIUnvCVUgza22zlNBRDUm9gf2AxfHDb/kGDID168PEY3XyROCcS1CcqqFT6nNgM9su6RvAo0AJMMPMFkiaHD0/DfgRcIekfxGqki4zs1X1OV9LlT5BTWlphg02bAg3TwTOuYRkG2KiW/RwXX0PbmaPAI/UWjct7fH7wEn1PX5rkD5BzfBMFXArV4Z7TwTOuYRkKxHMIzTu1tXoOziRiIpMeokgI+9M5pxLWLYOZYOaMpBi1bt3mHOmzktIPRE45xIWp7HYJaikJPQw9hKBc65QPBE0A1lnKkslgp49mywe51xx8UTQDGTtXbxiRZjcuPYE9s4510jiXDWUkZmtafxwik+fPmFyGth1KKHevWH5crwPgXMucXGvGhoAVEePSwkdwbwxuRGkkkCd61eu9ETgnEtUnVVDZjbIzAYTOoR9zsx6mFl34LPAH5sqwKLnJQLnXMLitBEcFnUMAyCam+CY5EJyu/BE4JxLWJwhJlZJ+j5wJ6Gq6MvA6kSjcoGZJwLnXOLilAjOBnoCfwIeBHpF61zSampg+3ZPBM65ROUsEURXB13SBLEUpd69MzcY9+6NdyZzzjWJnCUCSUMlTZc0S9ITqVtTBFcMli8PNUDdusHFF4fHZmmXjoInAudcouK0EdwHTANuBT5ONpziVVYG1dW1VnoicM41gTiJYLuZ/SrxSIpcaWloEtiFDy/hnGsCcRqLH5J0kaS9JXVL3RKPrMhkLRH06NHk8TjnikecEsG50f1309b5fASNrKwMqqpqrVyxIjQetGtXkJicc8UhzlVDPpREE6izRODtA865hMUpESDpIOBAoGNqnZn9LqmgilGdbQSeCJxzCcuZCCRdBYwlJIJHCJPZ/wPwRNCIyspgyxbYtAk6dYpWrlwJw4YVNC7nXOsXp7F4AnA8sNzMzgdGAh0SjaoIlZWF+12qh7xE4JxrAnESwSYz2wFsl9QVWIE3FDe63RLB9u2werUnAudc4uK0EcyVVAr8hjBHwXrghSSDKkalpeF+ZzvBqlXh3hOBcy5hca4auih6OE3S34CuZvZqsmEVn91KBN6r2DnXRGJdNZRiZksSiqPoeSJwzhWKT17fTHgicM4ViieCZmKvvcL9zjYCH2fIOddEPBE0E23bwp571ioRtG37SSuyc84lJO9EIOmN6PaNJAIqZjuHmaishBtuCJeQDh4clp1zLiF5NRYDmNkBkroDRyQQT1ErLYXh/6qE+ybBxo1h5bvvwqRJ4fHEiQWLzTnXesWZoWwPSW2ix0MljQc+MrO/xNh3nKSFkt6RdHkd24yVNF/SAkmz834FrUhZGUxcMOWTJJCycSNMmVKYoJxzrV6cqqGngY6S+gKPA+cDd+TaSVIJcDNhbKIDgbMlHVhrm1LgFmC8mQ0Dzsgn+NamrAx6bVma+cmldax3zrkGipMIZGYbgdOBG83sC4Qv9lxGA++Y2WIz2wrcDXy+1jZfAv5oZksBzGxF/NBbn7IyeL9kQOYnB9Sx3jnnGihWIpA0BpgIpKqD4rQt9AWWpS1XRevSDQXKJD0laZ6kc+oIYJKkuZLmrly5MsapW6bSUriyZCp07rzrE507w9SpBYnJOdf6xUkElwJXAH8yswWSBgNPxthPGdZZreW2wKHAacDJwJWShu62k9l0MxtlZqN6tuLr6svK4I6tE9n+s19+snLgQJg+3RuKnXOJiTPW0GxgdtryYuCbMY5dBfRPW+4HvJ9hm1VmtgHYIOlpwjDXb8U4fquT6l289rAT6A5w221wwQWFDMk5VwTqTASSHmL3X/A7mdn4HMd+ERgiaRDwHnAWoU0g3Z+BmyS1BdoDhwO/iBF3q5RKBOuXVYdEkFrhnHMJylYiuC66Px3oA9wZLZ8NLMl1YDPbHnU6exQoAWZEVUuTo+enmdkb0YimrwI7gFvN7LV6vZJWINWJeOP7NeGBJwLXgmzbto2qqio2b95c6FCKWseOHenXrx/t2rWLvU+diSCqEkLSj8zsM2lPPRRV4eRkZo8QprdMXzet1vJPgZ/GjrgVS33vb/6getcVzrUAVVVV7LnnnpSXlyNlaiJ0STMzVq9eTVVVFYMGDYq9X5zG4p5RAzEAUVVP622xLaDU9/62FVEi8HGGXAuyefNmunfv7kmggCTRvXv3vEtlcS4D/RbwlKTF0XI58PX8wnNxpL73t6/0EoFrmTwJFF59/gY5SwRm9jdgCHBJdNvPzB7N+0wup9T3vlXXQJs2YThS51xepk6dyrBhwxgxYgQVFRU8//zzBYlj/vz5PPLIJzXjM2fO5JprrgHgvPPO4/77799tn6eeeorPfvazTRZjSs4SgaTOwLeBgWb2NUlDJO1nZg8nH15x6dgx3FRTHYoH/uvKtVJ9+sCHH+6+vndvWL68/sd97rnnePjhh3nppZfo0KEDq1atYuvWrfU/YAPMnz+fuXPncuqppwIwfvx4xo/PdbFlYcRpI7gd2AqMiZargB8nFlGRKyuDknXVXi3kWrVMSSDb+rg++OADevToQYcOHQDo0aMH++yzD+Xl5axatQqAuXPnMnbsWABmz55NRUUFFRUVHHzwwaxbtw6Aa6+9luHDhzNy5EguvzyMl7lo0SLGjRvHoYceytFHH82bb74JhF/3kydP5uijj2bo0KE8/PDDbN26lR/84Afcc889VFRUcM8993DHHXfwjW98Mnr/Y489tss+tW3YsIELLriAww47jIMPPpg///nPACxYsIDRo0dTUVHBiBEjePvttxv2phGvjWBfMztT0tkAZrZJXhGYmNJSaLemBvp7InAt16WXwvz59ds3+o7eTUUFXH999n1POukkrr76aoYOHcoJJ5zAmWeeyTHHHFPn9tdddx0333wzRx11FOvXr6djx4789a9/5cEHH+T555+nc+fOrFmzBoBJkyYxbdo0hgwZwvPPP89FF13EE088AcCSJUuYPXs2ixYt4thjj+Wdd97h6quvZu7cudx0000A3HHHHbucO9M+6aZOncpxxx3HjBkzqKmpYfTo0ZxwwglMmzaNSy65hIkTJ7J161Y+/vjj7G9KDHESwVZJnYg6l0naF9jS4DO7jMrKoON7XiJwrj66dOnCvHnzeOaZZ3jyySc588wzd9bLZ3LUUUfx7W9/m4kTJ3L66afTr18/HnvsMc4//3w6R2N+devWjfXr1zNnzhzOOOOTAZK3bPnka/CLX/wibdq0YciQIQwePHhnaSGbXPvMmjWLmTNnct11oUvX5s2bWbp0KWPGjGHq1KlUVVVx+umnM2TIkLzeo0ziJIKrgL8B/SVVAkcB5zX4zC6jsjLovKUaSvvn3ti5ZirXL/dsdQpPPdWwc5eUlDB27FjGjh3L8OHD+e1vf0vbtm3ZsWMHwC6XVl5++eWcdtppPPLIIxxxxBE89thjmNluV97s2LGD0tJS5tdRzKm9fZxKk1z7mBkPPPAA++233y7rDzjgAA4//HD+8pe/cPLJJ3Prrbdy3HHH5TxfNlnbCKIJacoIvYvPA+4CRpnZUw06q6tTWRl02eYlAufqY+HChbvUmc+fP5+BAwdSXl7OvHnzAHjggQd2Pr9o0SKGDx/OZZddxqhRo3jzzTc56aSTmDFjBhujCaLWrFlD165dGTRoEPfddx8QvqRfeeWVnce577772LFjB4sWLWLx4sXst99+7LnnnjvbHDLJtE+6k08+mRtvvBGzMNLPyy+/DMDixYsZPHgw3/zmNxk/fjyvvvpqQ94yIEciMLMdwDfMbLWZ/cXMHjazVQ0+q6tT6V7GnjtqPBG4Vq137/zWx7V+/XrOPfdcDjzwQEaMGMHrr7/OD3/4Q6666iouueQSjj76aEpKSnZuf/3113PQQQcxcuRIOnXqxCmnnMK4ceMYP348o0aNoqKiYmfVTGVlJbfddhsjR45k2LBhOxtvAfbbbz+OOeYYTjnlFKZNm0bHjh059thjef3113c2FteWaZ90V155Jdu2bWPEiBEcdNBBXHnllQDcc889HHTQQVRUVPDmm29yzjkZR+/Pi1LZps4NpCuBTcA9wIbUejNb0+Cz18OoUaNs7ty5hTh1k/jRFRu58po92PGTa2hz+WWFDse52N544w0OOOCAQofR5M477zw++9nPMmHChEKHslOmv4WkeWY2KtP2cdoIUuMgX5y2zoDBGbZ1DdS7fehVvLlDKZ1zbOucc40hznwE8Ucucg3Wq11IBOvblXkicK4FqH1ZaEuUs0OZpM6Svi9perQ8RFLT94EuEt3a1ACwto23ETjnmkY+PYuPjJa9Z3GCuimUCGrkicA51zTiJIJ9zexaYBuEnsVkno/YNYKuO0IiWP1xaWEDcc4VjTiJwHsWN6E9t4VEsHK7lwicc00jTiL4Ibv2LH4c+F6SQRWzzltrAPhwS2lB43CuJerSpUvObZ555hmGDRtGRUUFmzZtyuv4Dz74IK+//vrO5R/84Ac89thjecfZ3MSZj2AW3rO4ybTfUM1aulL9UUnujZ1rySorobw8zL1RXh6Wm+S0lXznO99h/vz5dOrUKa99ayeCq6++mhNOOKGxQ2xyca4amgmcBDzlPYuTp5pqPmpTSnV1oSNxLkGVlTBpErz7LpiF+0mTGi0ZPPXUU4wdO5YJEyaw//77M3HiRMyMW2+9lXvvvZerr76aiRMnAvDTn/6Uww47jBEjRnDVVVftPMbvfvc7RowYwciRI/nP//xP5syZw8yZM/nud79LRUUFixYt2mWCmccff5yDDz6Y4cOHc8EFF+wclC7fIbALIU6Hsp8BZwLXSHqB0MP4YTPLb1JMF09NDevblnkicC1brnGo//lP2FKrqXHjRvjKV+A3v8m8T5xxqNO8/PLLLFiwgH322YejjjqKZ599lq9+9av84x//2NkTeNasWbz99tu88MILmBnjx4/n6aefpnv37kydOpVnn32WHj16sGbNGrp168b48eMz9iLevHkz5513Ho8//jhDhw7lnHPO4Ve/+hWXXnppnfFlGgK7UOJUDc02s4sIPYmnA18EViQdWNGqrmZDB08ErpWrnQRyra+H0aNH069fP9q0aUNFRQVLlizZbZtZs2Yxa9YsDj74YA455BDefPNN3n77bZ544gkmTJhAjx49gDAUdTYLFy5k0KBBDB06FIBzzz2Xp59+Ous+qSGwb7jhBmpqamjbNs7v8mTEOnN01dDnCCWDQ4DfJhlUUauuZnOnIdTUFDoQ5xog1y/38vJQHVTbwIENH4c6kpqlDMLQ1Nu3b99tGzPjiiuu4Otf//ou62+44Ya8JoHPNmZbPkNg77///rHP2ZjitBHcA7wBHAfcTOhX8F9JB1a0qqvZtoe3EbhWbupU6FxrEJXOncP6JnTyySczY8YM1q9fD8B7773HihUrOP7447n33ntZvXo1wM5ZyuoaWnr//fdnyZIlO2cZ+/3vf79zZrR8hsAulLg9i/c1s8lm9kQ0NLVLSk0N2/f0qiHXyk2cCNOnhxKAFO6nTw/rm9BJJ53El770JcaMGcPw4cOZMGEC69atY9iwYUyZMoVjjjmGkSNH8u1vfxuAs846i5/+9KccfPDBLFq0aOdxOnbsyO23384ZZ5zB8OHDadOmDZMnTwbIawjsQokzDHU74ELgM9Gq2cA0M9uWcGwZtephqLdtg/btmXX0jzjtue+zdWv2mZyca06KdRjq5ijfYajjlAh+BRwK3BLdDonWucYWFQPUrYzt28NFFM45l7Q4jcWHmdnItOUnJL1S59au/qJEUNKjdOfiHnsUMB7nXFGIUyL4OBpfCABJg4GPkwupiEWXCrXrGcYZ8nYC51xTiFMi+C7wpKTFhFFHBwLnJxpVsYq++Tvu7YnAtUxmltdll67x5Wr3zSTODGWPSxoC7EdIBG+amY8+moTom79T35AIvC+Ba0k6duzI6tWr6d69uyeDAjEzVq9enXcv5ZyJQNLFQKWZvRotl0n6ipndEmPfccAvgRLgVjO7po7tDgP+CZxpZvfn8wJalSgRdOlbmr7oXIvQr18/qqqqWLlyZaFDKWodO3akX79+ee0Tp2roa2Z2c2rBzKolfY1wBVGdJJUQOqCdSJjV7EVJM83s9Qzb/R/waF6Rt0ZREaDrQK8aci1Pu3btGDTIpzhvieI0FrdRWjkv+uJuH2O/0cA7ZrbYzLYCdwOfz7DdfwEP4OMXQXU1G+lE931C1/hvfSv0I5CgT58Cx+aca7XiJIJHgXslHS/pOMKcBH+LsV9fYFnaclW0bidJfYEvANOyHUjSJElzJc1t1cXO6mpqKM341IcfNm0ozrniEScRXAY8QehdfDHxZyjL1FpUuzn7euAyM8t6OaqZTTezUWY2qmfPnjFO3ULV1FCNT1HpnGtaca4a2iHpDuAJM1uYx7GrgP5py/2A92ttMwq4O6p56gGcKmm7mT2Yx3laj+pqTwTOuSYXZ/TR8cB8ouogSRXRrGW5vAgMkTRIUnvgLGCX/cxskJmVm1k5cD9wUdEmAfBE4JwriDhVQ1cRGn5rAMxsPlCeaycz2w58g9DG8AZwr5ktkDRZ0uR6xtu6ZWkjcM65pMS5fHS7ma2tTwcRM3sEeKTWuowNw2Z2Xt4naG1qatjSqQw27f5U795NH45zrjjEKRG8JulLQImkIZJuBOYkHFfx2bED1q7lq98twwyuuy6sXrs2zO29fHlhw3POtV5xEsF/AcOALYRLR9cClyQZVFFKfeOXhTaC1MVRK7x3hXMuYXEmr99oZlPM7LBoUoM7gZuSD63IpLoRl5YC0KtXWPRE4JxLWp2JQNIISbMkvSbpR5J6S3oAeAx4va79XD2lRpiLSgSeCJxzTSVbieA3wB+A/wesAl4CFgOfMrNfNEFsxSVVIvBE4JxrYtmuGupgZndEjxdK+g5wea5ewK6eaiUCbyNwzjWVbImgo6SD+WSoiPXAiNQAdGb2UtLBFZVabQQdOsBee3kicM4lL1si+AD4edry8rRlA45LKqiiVKuNAEL1UGseY8851zzUmQjM7NimDKToVVdD27a7zFbfq5eXCJxzyYvTj8A1herqUBpI68HticA51xQ8ETQXNTU72wdSPBE455qCJ4LmIlUiSNOrF6xaBR/7dVrOuQTlnQgk7S2pQxLBFLU6EsGOHbBmTYFics4VhfqUCH4PvCnpusYOpqjVkQjAq4ecc8nKOxGY2QnAYOD2xg+niGVoI/BOZc65phBnhrJ9U1VBksZK+iawl5ktSDy6YmHmJQLnXMHEKRE8AHws6VPAbcAgwhhErrHcfjts3w4/+QmUl0NlJeCJwDnXNOIkgh3RtJNfAK43s28BeycbVhGprIRvfOOT5XffhUmToLKSbt2gTRtPBM65ZMVJBNsknQ2cCzwcrWuXXEhFZsoU2FRrbsqNG2HKFEpKoEcPH2bCOZesOIngfGAMMNXM/i1pEGFyGtcYli7Nut47lTnnkhZnhrLXgcsI8xFgZv82s2uSDqxoDBiQdb0nAudc0uJcNfQ5YD7wt2i5QtLMhOMqHlOnQrtaNW2dO4f1eCJwziUvTtXQD4HRQA2Amc0nXDnkGsPEiXDEEVBSEgacGzgQpk8P6/FE4JxLXrb5CFK2m9lapY2KSZiPwDWWdu1g9GiYM2e3p3r1grVrYcuWMFmNc841tjglgtckfQkokTRE0o3A7t9Yrv6WLauzrSDVl8CvHHLOJSVOIvgvYBiwBbgL+Ai4NMGYiotZuEKof/+MT3unMudc0nJWDZnZRmBKdHONbeXKUO+To0TgicA5l5Q6E4Gkh8jSFmBm4xOJqNgsWxbuPRE45wokW4nAh5luCqkOZXVUDfkIpM65pGWbvH526rGk9sD+hBLCQjPb2gSxFYccJYI99wxXC3ljsXMuKXE6lJ0GLAJuAG4C3pF0SpyDSxonaaGkdyRdnuH5iZJejW5zJI3M9wW0eEuXQseO0L17xqcl70vgnEtWnH4EPwOONbN3IMxPAPwF+Gu2nSSVADcDJwJVwIuSZkZDVqT8GzjGzKqj5DIdODz/l9GCLV0aSgO79tPYhScC51yS4lw+uiKVBCKLgThfS6OBd8xscVSVdDfw+fQNzGyOmVVHi/8E+sU4buuybFmd7QMpngicc0nKdtXQ6dHDBZIeAe4ltBGcAbwY49h9gWVpy1Vk/7X/FeooZUiaBEwCGFDXIG0t1dKlcPLJWTfp1QsW+HxwzrmEZKsa+lza4w+BY6LHK4Gy3TffTaa6joyXo0o6lpAIPp3peTObTqg2YtSoUa1neItt2+CDD+oegTSSKhGYZa1Bcs65esl21dD5DTx2FZBe59EPeL/2RpJGALcCp5jZ6gaes2V5773w7R6jamjzZli/PlxF5JxzjSlb1dD3zOzaaGyh3X6Fm9k3cxz7RWBINJHNe8BZwJdqnWMA8EfgP83srXyDb/FyXDqakt6pzBOBc66xZasaeiO6n1ufA5vZdknfAB4FSoAZZrZA0uTo+WnAD4DuwC3R6KbbzWxUfc7XIuXoTJaSngj23TfhmJxzRSdb1dBD0cONZnZf+nOSzohzcDN7BHik1rppaY+/Cnw1drStTT0SgXPONbY4l49eEXOdy9eyZdCtG+yxR9bNPBE455KUrY3gFOBUoK+kG9Ke6gpsTzqwopDqTJaDjzfknEtStjaC9wntA+OBeWnr1wHfSjKoorFsGZSX59ysQwfo2tXHG3LOJSNbG8ErwCuS/kDoE+CDzjW2pUvh6KNjbeq9i51zSYkz1tCJwK8JA88JGCTp62aWdawhl8O6dVBTE6tqCDwROOeSEycR/Jx6DDrnckj1IchxxVCfPvDhh58sp3oW9+4Ny5cnFJtzrqgkOeicyyZ16WiOEkF6Eoiz3jnn8hWnRJBx0LnUoHRm9scE42u9YpYInHMuaXESQUd2H3SuG2FQOiMMEeHytXQptGkD++xT6Eicc0UuZyJohMHnXCZLl0LfvtA2Ti52zrnkxJmqcqikxyW9Fi2PkPT95ENr5WJMSOOcc00hTmPxbwhDSmwDMLNXCSOJuoaI2au4d+/81jvnXL7iJILOZvZCrXU+xERD7NgBVVWxSgTLl4cpC8zgllvCusWL/dJR51zjiZMIVkV9BwxA0gTgg0Sjau1WroQtW2J3Jks58shwP2dOAjE554pWnERwMaFn8f6S3gMuBS5MMqhWr56Xjh50EHTpAs89l0BMzrmiFeeqocXACZL2ANqY2brkw2rFKivh0kvD4wsvDPNPTpwYa9eSEjj8cC8ROOcaV5yrhv5XUqmZbTCzdZLKJP24KYJrdSorYdIkWLUqLH/wQViurIx9iCOPhFdeCfnDOecaQ5yqoVPMrCa1YGbVhHkKXL6mTIGNG3ddt3FjWB/TmDGhrfnFFxs5Nudc0YqTCEokdUgtSOoEdMiyvatLanyhuOszOOKIcO/VQ865xhKnW+udwOOSbidcOXQB8NtEo2qtBgyAd9/NvD6msjI48EBPBM65xpOzRGBm1wI/Bg4AhgE/ita5fE2dGlp803XuHNbnYcwY+Oc/QxWRc841VJzG4j2AWWb2HWA60EFSu8Qja40+97kwoUCXLuF+4ECYPj32VUMpRx4Ja9bAW28lFKdzrqjEaSN4GugoqS/wGHA+cEeSQbVaf/oTbN8Os2aFn/NLluSdBMA7ljnnGlecRCAz2wicDtxoZl8ADkw2rFbqzjth8OBPWnzraejQ0FbgHcucc40hViKQNAaYSJiiEuI1Mrt0H3wATzwRSgCp+SbraZ99oLoabr01HCp169OnkWJ1zhWVOIngUsLoo38yswWSBgNPJhpVa3T33aE6qB5VQbX59JXOucYkMyt0DHkZNWqUzZ07t9Bh5O/QQ8PP9kaIPVuBooX9OZ1zTUTSPDMblem5OFcNPSnpidq3xg+zlaqsDDORvfRSGD86j+EknHOuKcSp6/9O2uOOwP/D5yOIJzW2UGpYierqsAyNUkXknHONoV5VQ5Jmm9kxubdsfC2qaqi8PHNP4oEDw6Wj9RS3rbl3b5/AxjkXNLRqqFvarYekkwG/PqUulZUhAbRpkzkJQF5jC2USd5rKDz/0q4qcc7nFuWpoHjA3un8O+G/gK3EOLmmcpIWS3pF0eYbnJemG6PlXJR2ST/CxpX85l5eH5drrLroo+3LcfSZNCgkgW0krz5nJakufvjKfAl3txOA3v/mt5d4a9YedmSVyA0qARcBgoD3wCnBgrW1OBf4KCDgCeD7XcQ899FDLy513mnXunP69adaunVn79ruuy3Wrzz6Zbp07h5gaUUND8pvf/NYyb/l9TzDXLPP3ap0lAkmHSeqTtnyOpD9Hv+C7xcgxo4F3zGyxmW0F7gY+X2ubzwO/i+L8J1Aqae9YGSyuTHMAbNsGW7fmd5z67JNOqvfYQs45l6RsVUO/BrYCSPoMcA3wO2AtYfC5XPoCy9KWq6J1+W6DpEmS5kqau3LlyhinTtPA+vhGMXBgg8YWcs65JGVLBCVmtiZ6fCYw3cweMLMrgU/FOLYyrLN6bIOZTTezUWY2qmfPnjFOnaaB9fF5U62XVI9hpvMVt/HYOecyyZoIJKX6GRwPpHcii9P/oAron7bcD3i/Hts0zNSp4cs4Xbt20L59fseJs0/nzjB5cigBNGFVUO3GY08Mzrl8ZEsEdwGzJf0Z2AQ8AyDpU4TqoVxeBIZIGiSpPXAWMLPWNjOBc6Krh44A1prZB/m+iKwmTgxfxulfzrffDjNm7LruwguzL8fZZ/p0uOWWUAVUwKqg2onBk4NzrU9j/k9n7VAWfTnvTZiYZkO0bijQxcxeynlw6VTgesIVRDPMbKqkyQBmNk2SgJuAccBG4Hwzy9pbrEV1KHPOuWYiW4cyH3TOOeeKQIN6FjvnnGvdPBE451yR80TgnHNFzhOBc84VuRbXWCxpJfBuPXfvAaxqxHCS1pLibUmxQsuKtyXFCi0r3pYUKzQs3oFmlrFHbotLBA0haW5drebNUUuKtyXFCi0r3pYUK7SseFtSrJBcvF415JxzRc4TgXPOFbliSwRxRk1tTlpSvC0pVmhZ8bakWKFlxduSYoWE4i2qNgLnnHO7K7YSgXPOuVo8ETjnXJErmkQgaZykhZLekXR5oeOpTdIMSSskvZa2rpukv0t6O7ovK2SMKZL6S3pS0huSFki6JFrf7OKV1FHSC5JeiWL9n+Yaa4qkEkkvS3o4Wm7OsS6R9C9J8yXNjdY153hLJd0v6c3o8zumOcYrab/oPU3dPpJ0aVKxFkUikFQC3AycAhwInC3pwMJGtZs7CMNxp7sceNzMhgCPR8vNwXbgv83sAOAI4OLo/WyO8W4BjjOzkUAFMC4aXr05xppyCfBG2nJzjhXgWDOrSLu+vTnH+0vgb2a2PzCS8D43u3jNbGH0nlYAhxKG6f8TScVa16z2rekGjAEeTVu+Arii0HFliLMceC1teSGwd/R4b2BhoWOsI+4/Ayc293iBzsBLwOHNNVbCLH2PA8cBDzf3zwGwBOhRa12zjBfoCvyb6CKZ5h5vWnwnAc8mGWtRlAiAvsCytOWqaF1z19uiGdui+14Fjmc3ksqBg4HnaabxRlUt84EVwN/NrNnGSpjI6XvAjrR1zTVWCHOMz5I0T9KkaF1zjXcwsBK4Pap6u1XSHjTfeFPOIswYCQnFWiyJQBnW+XWzDSSpC/AAcKmZfVToeOpiZh9bKGL3A0ZLOqjAIWUk6bPACjObV+hY8nCUmR1CqHa9WNJnCh1QFm2BQ4BfmdnBwAaaQTVQNtE0v+OB+5I8T7Ekgiqgf9pyP+D9AsWSjw8l7Q0Q3a8ocDw7SWpHSAKVZvbHaHWzjRfAzGqApwhtMc0x1qOA8ZKWAHcDx0m6k+YZKwBm9n50v4JQhz2a5htvFVAVlQgB7ickhuYaL4QE+5KZfRgtJxJrsSSCF4EhkgZFGfYsYGaBY4pjJnBu9PhcQl18wUVzTd8GvGFmP097qtnFK6mnpNLocSfgBOBNmmGsZnaFmfUzs3LCZ/QJM/syzTBWAEl7SNoz9ZhQl/0azTReM1sOLJO0X7TqeOB1mmm8kbP5pFoIkoq10A0hTdjgcirwFrAImFLoeDLEdxfwAbCN8MvlK0B3QsPh29F9t0LHGcX6aULV2qvA/Oh2anOMFxgBvBzF+hrwg2h9s4u1Vtxj+aSxuFnGSqhzfyW6LUj9XzXXeKPYKoC50efhQaCsucZLuLhhNbBX2rpEYvUhJpxzrsgVS9WQc865OngicM65IueJwDnnipwnAuecK3KeCJxzrsh5InCuDpK6p43+uFzSe9Hj9ZJuKXR8zjUWv3zUuRgk/RBYb2bXFToW5xqblwicy5OksWlzBfxQ0m8lzYrG5j9d0rXRGP1/i4biQNKhkmZHg7M9mhomwLnmwBOBcw23L3Aa8HngTuBJMxsObAJOi5LBjcAEMzsUmAFMLVSwztXWttABONcK/NXMtkn6F1AC/C1a/y/CHBP7AQcBfw/DNFFCGE7EuWbBE4FzDbcFwMx2SNpmnzS87SD8jwlYYGZjChWgc9l41ZBzyVsI9JQ0BsIQ3pKGFTgm53byROBcwsxsKzAB+D9JrxBGaz2yoEE5l8YvH3XOuSLnJQLnnCtyngicc67IeSJwzrki54nAOeeKnCcC55wrcp4InHOuyHkicM65Ivf/A0y3CSbpSADyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.integrate as spi\n",
    "import numpy as np\n",
    "import pylab as pl\n",
    "\n",
    "beta=1.4247\n",
    "TS=1.0\n",
    "ND=70.0\n",
    "S0=1-1e-6\n",
    "I0=1e-6\n",
    "INPUT = (S0, I0)\n",
    "\n",
    "def diff_eqs(INP,t):\n",
    "\t'''The main set of equations'''\n",
    "\tY=np.zeros((2))\n",
    "\tV = INP\n",
    "\tY[0] = - beta * V[0] * V[1]\n",
    "\tY[1] = beta * V[0] * V[1] \n",
    "\treturn Y   # For odeint\n",
    "\n",
    "t_start = 0.0; t_end = ND; t_inc = TS\n",
    "t_range = np.arange(t_start, t_end+t_inc, t_inc)\n",
    "RES = spi.odeint(diff_eqs,INPUT,t_range)\n",
    "\n",
    "print(RES[:3])\n",
    "\n",
    "#Ploting\n",
    "pl.plot(RES[:,0], '-bs', label='Susceptibles')  # I change -g to g--  # RES[:,0], '-g',\n",
    "pl.plot(RES[:,1], '-ro', label='Infectious')\n",
    "pl.legend(loc=0)\n",
    "pl.title('SIR epidemic without births or deaths')\n",
    "pl.xlabel('Time')\n",
    "pl.ylabel('Susceptibles, Recovereds, and Infectious')\n",
    "#pl.savefig('2.1-SIR-high.png', dpi=900) # This does, too\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-23T14:25:58.338419Z",
     "start_time": "2021-09-23T14:25:58.313168Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Susceptibles</th>\n",
       "      <th>Infectious</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>9.999990e-01</td>\n",
       "      <td>0.000001</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>9.999958e-01</td>\n",
       "      <td>0.000004</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>9.999825e-01</td>\n",
       "      <td>0.000018</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>9.999271e-01</td>\n",
       "      <td>0.000073</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>9.996971e-01</td>\n",
       "      <td>0.000303</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>66</th>\n",
       "      <td>-6.352040e-13</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>67</th>\n",
       "      <td>-5.844091e-13</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>68</th>\n",
       "      <td>-5.336141e-13</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>69</th>\n",
       "      <td>-4.828191e-13</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>70</th>\n",
       "      <td>-4.320242e-13</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>71 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "    Susceptibles  Infectious\n",
       "0   9.999990e-01    0.000001\n",
       "1   9.999958e-01    0.000004\n",
       "2   9.999825e-01    0.000018\n",
       "3   9.999271e-01    0.000073\n",
       "4   9.996971e-01    0.000303\n",
       "..           ...         ...\n",
       "66 -6.352040e-13    1.000000\n",
       "67 -5.844091e-13    1.000000\n",
       "68 -5.336141e-13    1.000000\n",
       "69 -4.828191e-13    1.000000\n",
       "70 -4.320242e-13    1.000000\n",
       "\n",
       "[71 rows x 2 columns]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df = pd.DataFrame(RES, columns = ['Susceptibles', 'Infectious'])\n",
    "df.to_excel('Fig2.3-SI.xlsx', index = True)\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-23T12:42:26.681600Z",
     "start_time": "2021-09-23T12:42:25.438397Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[9.99999000e-01 1.00000000e-06 0.00000000e+00]\n",
      " [9.99996055e-01 3.64966741e-06 2.95304063e-07]\n",
      " [9.99985485e-01 1.31593559e-05 1.35516202e-06]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABPJklEQVR4nO2dd3yUVfa4n5OEDoIUQaUroNJCsbCIIKBiw9WNisuuiosssmvdYmFdVGR/fl13VWzYAFdib4uKDZBqWxBUqkBogdBbqCGZ8/vjfWeYTGYm7ySZZJI5D5/3M3Pve+/7npkM97z3nHvPEVXFMAzDSF5SKloAwzAMo2IxRWAYhpHkmCIwDMNIckwRGIZhJDmmCAzDMJIcUwSGYRhJjikCo1hE5BMRuSHCudYioiKSFof7ThCR+8v6ukHX3y8ibaOcXyciA+N1/2iIyCwRGR7hXEtX9tQYr6kicmrZSFh2RPusiXztqoQpggRCRM4Vka9EZK+I7BKR+SJypnvuRhGZF9R2nYgccgeELSIyWUTqxkMuVb1YVV+Jx7WLue9IVR0bx+vXVdUsAPf7ezhe9wom9G8ZK6q6wZW9IMo9km4AFJEHRGRKRctRGTFFkCCIyHHAR8BTQEPgZOBB4EiUbperal0gHegG3BtnMY0ERxwS4v91PGaJRnxIiB+MAUB7AFV9XVULVPWQqn6uqj8W11FVtwCf4SiEsIhIfRF5WURyRGSTiDzsNy24T6jzReQpdzayQkQGBPUNPF2KSKqIPCYiO0QkC7i0BPd5XET2iEiWiPzCrd8oItuCTVChT+kicoWILBaRfSKyRkQGhfmcw0Tkw6DyahF5K6i8UUTS3fcqIqeKyAhgKPBXd4b1YdAl00XkR/d7eVNEagZd62b3+rtEZKqInOTWFzGX+b9DETkdmAD0cu+1J9LfDDhFRL5z7/1fEWkY7vrutceJyHzgIPAq0Ad42r3H00HXHCgiq0Rkt4g8IyLiXuNUEZnt3muHiLwZSSgRGSwiS92/4Sz3M/nPrRORu0XkR+BAOGUgIhe4v7G9rmwScv4mEVnuyviZiLQKOvek+zfcJyILRaSPWz8IuA+41v3MPwRdspX7u8sVkc9FpLHbp6aITBGRne5n+Z+INI3y96i6qKodCXAAxwE7gVeAi4HjQ87fCMwLKq8DBrrvmwM/AU9Guf4HwPNAHeAE4Dvg90HXzgfuBKoB1wJ7gYbu+VnAcPf9SGAF0AJn5vIloEBaDPcZBqQCDwMbgGeAGsCFQC5Q120/GXjYfX+WK9MFOA8wJwOnhfmcbYE9bpsTgfXApqBzu4EUt6zAqaH3CvmOvwNOcj/rcmCke64/sAPo7sr+FDDHPdc6+DsJ8x0W+ltG+HvNAjYBndzv8l1gSrjru203AB2BNPdvGLhf0DUVZ9bZAGgJbAcGuedeB0a731tN4NwIcrUHDrh/h2rAX4HVQPWg72wxzu+jVpj+jYF9QIbb/073N+H/bn7pXu9097P8DfgqqP9vgEbuuT8BW4Ca7rkH/N9RyPe4xpW7llt+xD33e+BDoDbO77EHcFxFjwUVcdiMIEFQ1X3AuTj/WV8EtrtPmdGeUD4QkVxgI7ANGBOukXuNi4E7VPWAqm4DHgeGBDXbBjyhqkdV9U1gJSFP+y7XuO02quou4P/FeJ+1qjpJHfv2mzgDxkOqekRVPwfygHAOzd8BE1X1C1X1qeomVV0R2kgdm38uzuyoL85MaZOInOaW56qqL9z3FIHxqrrZ/awfcmzWNdSV53tVPYJjluslIq1juHZxvKqqS1T1AHA/cI1EdhBPVtWlqpqvqkejXPMRVd2jqhtwlHi6W38UaAWcpKqHVTWSD+Na4GP373AUeAxngP1FUJvx7u/jUJj+lwDLVPUdt/8TOIO5n98D/09Vl6tqPvAPnFlZKwBVnaKqO93P+S8cJdwhyucFmKSqP7vyvBXymRvhPAwUqOpC9/9h0mGKIIFwf/w3qmpznCfBk3D+o0Til6paD+gHnIbztBWOVjhPXznuFHgPzlP7CUFtNqn7mOSy3r1/KCfhKJ7gdrHcZ2vQ+0MAqhpaF87p3QLnyc4Ls3G+k/Pc97NwlEBftxwLwYPUwSDZTiLos6vqfpwZ3ckxXj8aod9zNSL/jTdGqA8l0uf5K46J5jvX7HNThP6hn9vn3jv4c0eTpdDvx/3NBbdvBTwZ9PvZ5cp1MoCI/Mk1G+11z9cn8nfiJ9JnfhXnQeENEdksIo+KSLVirlUlMUWQoLhPu5NxFEJxbWe7bR+L0GQjjtO5sao2cI/jVLVjUJuT/fZil5bA5jDXysEZlIPbxXKfkrIROMVjW78i6OO+n03xiiDWMLybcQYtAESkDs7T5SYc0wk4Jgc/zUpwr9Dv+SiOOSocodeM6fOo6hZVvVlVT8J5Kn9Wwi81Df3c4sq5yeO9C/1+gvr72YhjSmwQdNRS1a9cf8DdOLPS41W1AY650P+7jfUzH1XVB1X1DJwZzWXA9bFco6pgiiBBEJHT3Ked5m65BXAd8I3HSzwBXOB3hAajqjnA58C/ROQ4EUkRkVNEpG9QsxOA20SkmohcjWOjnRbmPm+57ZqLyPHAPTHep6S8DAwTkQHudU92zT3hmA2cj2OjzgbmAoNwBupFEfpsxfEheOU1V550EamBY8L4VlXXqep2nIHxN+I412+isBLbCjQXkerF3OM3InKGiNQGHgLe0ShLRkvzeUTkav9vD8ePokC4e70FXOr+Harh2OmPAF95vNXHQEcRucp1JN9GYSU5AbhXRDq6ctV3f48A9XD8CduBNBH5O45vzc9WoLV4XDUlIueLSGfX3LYPR9F6/X6rFKYIEodc4GzgWxE5gKMAluD8RysWd/D5D44tORzXA9WBZTj/0d/Bcab6+RZoh/PEOQ7IUNWdYa7zIs50+gfge+C9GO9TIlT1Oxwn8+M4T4GzCXoyDWn7M7AfRwH4/S9ZwPwoA+nLwBmuSeIDD/LMwPmu38V5yj2Fwr6Qm4G/4JiLOlJ4oJwJLAW2iEikJ3xwTBeTcR2iOIOmV54EMtyVN+M9tD8T57e3H5gK3K6qa0MbqepKHIftUzi/lctxljHneRFKVXcAVwOP4Hw37YD5QeffB/4Px1yzD+f/wMXu6c+AT4CfccxThylsVnrbfd0pIt97EKcZzu9zH85CgNlAUu5DkMJmYSMZEZEbcVZtnFvRshiGUf7YjMAwDCPJMUVgGIaR5JhpyDAMI8mxGYFhGEaSU+mCQjVu3Fhbt25d0WIYhmFUKhYuXLhDVZuEO1fpFEHr1q1ZsGBBRYthGIZRqRCR9ZHOmWnIMAwjyTFFYBiGkeSYIjAMw0hyTBEYhmEkOaYIDMMwkpy4KQIRmShO6sElEc6LiIwXJ9XfjyLSPR5yNGsGIuV3NGtWvEyGYcSfnNwc+k7uy5b9WzyVK1ufsiSeM4LJOKF/I3ExTuTBdsAI4Ll4CLF1a/Ftyvp+phyMqkA8Brh4XTdcn7FzxjJvwzzGzh7rqVzZ+pQpxeWyxAmvW8N93w8nFG4DL3kwcXKrLolw7nnguqDySuDE4q7Zo0cPjQVIvKNp05g+gmHExOZ9m/W8SedpTm5OxLriyqqqt3x0i6Y8mKKjPhrlqVzSNiM/HKkpD6boiKkjdN/hfXrTBzdpygMpesP7N2j23mxdt3udDn13qMoDote9c50u2bpEr37rapUHRDPezNDvsr/TK9+4UuUB0V++/kudvW62vrP0Ha0+trryAFp9bHV9+tunC5Uf/+rxQuUXFrygE/43Qas9VM2pe6i6Pjz74UC52kPV9MFZDxYqPzb/MX103qOF6u6dfm+h8t1f3F2oPHb2WH1o1kOa9lBaoO5Pn/6pUPnOT+4sVL5/xv36txl/C9TVerhWob+TV4AFGmFcLTbWkIgsBnq6g/pnOLHKO6jqJcUpGTd/60eqWiTLloh8hJM/dZ5bngHcrapFdouJyAicWQMtW7bssX59xH0RYWTw3LRcsRBPhhdycnMY8u4Q3sx4k2Z1m4WtCy2P+ngUzy98npE9RvLMpc8AFKnzl4d3H87D5z/MnZ/dyWs/vcZVp13FH876Axv3bWT41OEc9R2lWko1bj3rVsZ/N558Xz5pKWlknJ7B28vepkALSJVUBrYdyJGCI8xZPwef+hCEUxueyupdq1EUQWhSuwn5ms+uQ7sq8iut9FRPrc7wbsMDf1uviMhCVe0Z9pwHRfC9qnYXkb8Ah1X1KRFZpKrdPNy4NZEVwcc4SaqDFcFfVXVhtGv27NlTY9lZbIrASGRKO6g/Puhxhk8dzpQfpzC4w2Aua3cZo6aN4qjvKGkpafy606/Ze2QvH/78YWCAblqnKVsOlM7OnEIKPnyBcqNajUiRFHYc3BEY+OtUr8PBvIP48JFCCh0adSAlJYXlO5bjUx+pksrxtY5n96HdAYXSuHZjdhzcESif2/JcRIS56+dSoAWkpaTRrE4zthzYElBKJ9U9ic37NwfK/dv0Z9a6WeQVeMqVA0D1lOogxNSnRmoNAI4UHCn3PrXSapF1e1bg4cALpVUE3+KkQRyNk4lorYgsCTe4h+nbmsiK4Hlglqq+7pZXAv3USXcYEVMERmUl3NN9pCf1YenDGNp5KIMyB5FXkEdaShq/6fwbth7YyqerP0U9puetmVaTGqk12HdkH4qSQgr1atQjNy83MBg3qdOE7Qe2BwbaAW0GMGvdrJgGq5qpNQE4XHC4VG3Kok+qpJKakhrToJ5CCgj41Fd84wToU5JZQTRF4MVZPAzoBYxzlUAbyiad21Tgenf10DnA3uKUQElo2rSsr1g2mDO56hPqwAx29h06eoiZWTN5+fuX8amP5xc+T5+JfZiwYAI+9fHyopfp/5/+gcEs35fPa0te45vsYymsUySFJrWbkCqpgDMApoSm61U4kn8koDh8+Nh7ZG9gUCnQArbs30KBm8Ez35fP9KzpnhWNn7yCPPJ80QdeL23Kok+BFsSkBMD5XmIZnCu6T15BHl9le00TXTxxy0cgIq/jOJcb4ySVHgNUA1DVCSIiwNM4K4sOAsPC+QdCiXVGUJ40a1byVUo2Q6j8hDXrLHieKzpcwTnNz+G+mfdRoAUIEnagrZFag7yCvMCTe+hTYEmejkvyBFoVSG+WzqLfL6poMRKK0pqG1kLRX62qti0b8WIjkRVBKLEoBlMElYtwZp4b3r+BV398lZ4n9aRhrYZ8tuazsH1TSKFz084s276Mo76jnu9ZnoO6DaRVj2iKwEsY6uCONYGrgYZlIVhVZ0uIPy5R/RVG8YQO/A/OfpC56+dy9VtX07x+c+aun8um3E0A/G/z/ziu+nEBh6rfdOM3v/jwsWTbElJTUmOSwYcvzCNZUWwQN2KlWEWgqjtDqp4QkXnA3+MjkmFUPKED/0OzH2Lu+rlkvJVB7Wq1+SLrCwDmbZxH011NqVO9DqmSSoEWUC2lGgeOHgisqvErgGAKtICCgqL1odigbpQHxSqCkNAPKTgzhHpxk8gwEgC/Y/em/97E8TWP57UlrwEwf+N86lavS4qk4FMf1VOqc9EpF/HWsrcCA75Xc48N8kai4MU09K+g9/nAOuCauEhTxWnaNLzPIFFXNiULwU//tavVZvLiyTy/8Hl86uOT1Z+QJmkBM0+apHHo6KGAnT7Pl0fmT5mezDw28BuJihfT0PnlIUgyEOwzePddyMiAWbOgb98KEykpCTX7jJk1hrnr59L75d5s3r+Zw/nHVuWkSiqKBsw8+Zpf5HqRzDw28BuVBS+mofo4Sz/Pc6tmAw+p6t54ClbVufBCqF4dpk41RVDe+M0+d356J/Vq1OPF718EIGtPFkM6DuG9Fe8F1qGHs++HwwZ9ozLjZfnou8AS4BW36rdAV1W9Ks6yhaUyLR8tjksugZ9/hlWrbEVRPAmeAeTl53HqU6cG7Pji/vPho3pqddo1bMeqXauK3ZBkA79R2Sjt8tFTVPVXQeUH3UB0RikZPBhuuQWWL4czzqhoaaou/hlAxlsZ/LD1h4ASSCEFEQk89ecV5LFs+7Kwm71s4DeqMl4UwSEROTcoOFxv4FB8xUoOLrvMUQRTp5oiiBfLty/nxe9fxKc+5m+cj3Bs6hVuXX611GoliuxoGJUZL7GGbgGeEZF1IrIeJyzEyPiKlRw0bw49esB//1vRklQd/PF9snZn8ci8R0h/Pp18n+PgFQQpxgZX1jFcDKMy4GXV0GKgq4gc55b3xVuoZOKKK2DMGGdFkQWfKz3+Hb+dn+3MwfyDhYKwKYEkSIUws4+R7EScEYjIb9zXu0TkLmA4MDyobJQBgwc7cYY+/riiJamcBEf4fG/Ze7yw8AUU5XDBYQadMoi0lMLPOtVTqzOq5yh0jAYOUwJGshPNNFTHfa0X5qgbZ7mShosucl6HD7fQ1CVh7JyxzFs/j36T+/Grt38VcPSmpaSxYPOCIqt/zPRjGEWJaBpS1efdt9NVdX7wOddhbJQBkaKTljScdTKxed9mXvr+JXz4WLlzZSDWDzgD/oGjB8j5U05MWZwMIxnx4ix+ymOdYcSVYDPQtgPb+MXEXxTaDxC67LNACxg7e2xFiGoYlYqIMwIR6QX8AmgS4hM4Dogtfq5hlAH+/QDDpw7nm+xv2HnoWGDccI5gMwMZhjeirRqqjuMLSKNwtNF9QEY8hTKMUHJyc5i0eBI+9fHxqo9pUKMB1VKqFYr0WZI8roZhRPcRzAZmi8hkVV1fjjIZRhHunn43R/KdZOr+5OSh4Z5tBmAYJcPLzuKXRORqVd0DICLHA2+o6kVxlSxJsNDUkfHHCBrebTiv/vhqoL5ACzh49KA5gg2jjPDiLG7sVwIAqrobOCFuEiUZW7Y4+wiuvx5atHDeqxZNc5mMjJ0zlrnr53L9B9cXOWeOYMMoO7woAp+ItPQXRKQVnjKnGrHQsiVs2gT5RcPdJyWb923mxe9fDBsADswMZBhliRfT0GhgnojMdsvnASPiJ1Jy0rIl+HywebPzPpkp8BVwwasXBGIEmRPYMOJLsTMCVf0U6A68CbwF9FDVz+ItWLLhH/zXJ7FbPic3h/Mmncdlr13Gsh3LAvV5BXlMWjyJLfvNXmYY8aBYRSBOuMZBQHdV/RCoLSJnxV2yJKNVK+d1w4aKlaMiGTNrDHM3zOXTNZ+SKoW3qphPwDDihxcfwbNAL+A6t5wL2By9jGnRwnlNVkWQtTuLlxe9HCiHpog0n4BhxA8vPoKzVbW7iCwCZ9WQiFSPs1xJR5060KhRciqCowVH6f9Kf3zqJIg3n4BhlC9eZgRHRSQVd6WQiDQBfHGVKklp2TL5FIFPfVz37nWs33vMOWI+AcMoX7wogvHA+8AJIjIOmAf8I65SJSktWyaPszgnN4e+k/oy8sORvLv8XfMJGEYFEi3oXBtVXauqmSKyEBgACPBLVV1ebhImEa1awYwZzoayYjIqVnrGzhnL3A1zmbNhDo1rN2bHwR2FzptPwDDKj2g+gneAHiIyQ1UHACvKSaakpWVL2L8f9u6FBg0qWpr4kZObw0vfv4SipEoqP4z8gZPqnVTRYhlG0hJNEaSIyBigfbjUlKr67/iJlZz49xJs2FC1FcGtn9waCBiXmpLKuDnjzDFsGBVINB/BEOAwx8JQhx7FIiKDRGSliKwWkXvCnK8vIh+KyA8islREhsX+EaoOwYqgqvLT1p94d/m7gbI5hg2j4okWhnol8H8i8qOqfhLrhd2VRs8AFwDZwP9EZKqqLgtq9gdgmape7q5GWikimaqaF+aSVZ6qrAhycnO49p1r+Xnnz0XO+R3DNiswjIrByz6CmSLya6B1cHtVfaiYfmcBq1U1C0BE3gCuAIIVgQL13N3LdYFdQNKGXWvaFKpXr5orh/zO4XCYY9gwKhYviuC/wF5gIXAkhmufDGwMKmcDZ4e0eRqYCmzGMTddq6pF9iiIyAjcQHctq3BEtpQUZ4dxVZsR+J3DAGmSxsa7NloeAcNIILwoguaqOqgE1w63ADI0pvBFwGKgP3AK8IWIzFXVfYU6qb4AvADQs2fPKh0CuypuKvvLF38JOIdTUlLMDGQYCYaXDWVfiUjnElw7G2gRVG6O8+QfzDDgPXVYDawFTivBvaoMVU0RZO/N5vWfXg+UzTlsGImHF0VwLrDQXf3zo4j8JCI/euj3P6CdiLRxYxMNwTEDBbMBZ6MaItIU6ABkeRe/6tGypZOT4OjR4ttWBjLezsAXEpHEdg0bRmLhxTR0cUkurKr5IvJH4DMgFZioqktFZKR7fgIwFpgsIj/hmJLuVtUdES+aBLRq5SSo2bQJWreuaGlKTk5uDpe9fhnf53xf5Jw5hw0jsYgWYqKh+za3pBdX1WnAtJC6CUHvNwMXlvT6VZHgJaSVWRGMmTWG73O+p3ZabTbcuYFGtRtVtEiGYUQg2oxgIY5zN5LTt21cJEpyqsJegpzcHCYumghAvuYHHMWGYSQmEX0EqtpGVdu6r6GHKYE4URUS1Nw67dZCiWXMH2AYiY0XZ7FRjtSuDY0bV15FkL03m/dWvBco2yohw0h8TBEkIK1aVV5FMPT9oWjIdhFbJWQYiY0pggSksiao2XFwB/M3zC9Sb6uEDCOx8bJqKCyquqvsxTHAUQRffFH5EtTcN+M+AJbcsoSOJ3SsYGkMw/CK11VDLYHd7vsGOBvB2sRbuGTFn6Bmzx44/viKlqZ4cnJzuOS1S1i8ZTF/6vUnUwKGUckodtUQzoawy1W1sao2Ai4D3ovUzyg9rVo5r5XFT/Dg7AdZvGUxtavVZkzfMRUtjmEYMeLFR3CmuzEMADc3Qd/4iZTcNGsGGRnO+/R0xzQk4tQnIoX2DBTkc+DogQqWyDCMWPGiCHaIyN9EpLWItBKR0cDOeAuWrGzdGlt9RTNm1phjG8bE9gwYRmXEiyK4DmgCvA98AJzg1hlJTk5uDpMWTwqUbc+AYVROilUEqrpLVW9X1W7ucbutGDIA7v/yfvJ9hRPK2Z4Bw6h8FBt9VETaA3+maKrK/vETy6gMTFs1rUid7RkwjMqHlzDUbwMTgJeAgmLaGknC7kO7OXj0IJe3v5yp14WmmTAMozLhRRHkq+pzcZfEAJwE9uEcw02blr8skcjJzeHsl85m75G9PNz/4YoWxzCMUuLFWfyhiIwSkRNFpKH/iLtkScqWLc6O4jvvhDp1nPeqTn2icO+Me9m4byOnNjyVLk27VLQ4hmGUEi8zghvc178E1Vk+gjjTqBEcOABHjkCNGhUtzTFycnOY8uMUADbu3ciW/VtoVjdBNzkYhuEJL6uGLB9BBdDITei1M8F2bIyeOTqQa0BRWyFkGFUAT9FHRaSTiFwjItf7j3gLluzErAgyM53clikpzmtmZpnLlJObw6s/vhoo274Bw6gaFKsIRGQM8JR7nA88CgyOs1xJT0yKIDMTRoxwYlerOq8jRpS5Mhgza4ztGzCMKoiXGUEGMADYoqrDgK5AAlmtqyaNGzuvnhTB6NFw8GDhuoMHnfoy5NPVnxaps30DhlH58eIsPqSqPhHJF5HjgG2YozjuFDsjyMx0BvoNG5xZQDjKMHypT33UqV6Hbs26sXDEQqQyJUowDCMqXhTBAhFpALyIk6NgP/BdPIUyjimCHTvCnPSbgkJnAaG0bFlm8kxbNY0VO1aQeVWmKQHDqGIUqwhUdZT7doKIfAocp6o/xlcso2ZNJ5F92BlBOFNQKLVqwbhxZSbPY189RovjWnD1GVeX2TUNw0gMYspZrKrrTAmUH40aRVAE0Uw+/qf1a66BoUPLRI5PVn3C7PWzuanbTVRLrVYm1zQMI3Gw5PUJTERFEMnk06oV+Hxw9tkwfz4UlE1oqFs/uRWA7H3ZZXI9wzASC1MECUxERTBuHKSFWPVq1z5mCrrrLli9Gj76qNQyLNy8kDW71wDw2k+v2Z4Bw6iCmCJIYBo3jqAIhg6FLl0cZSDizAReeOGYKeiqq5xZw+OPl1qGkR+NDLy3PQOGUTWJWRGIyHL3+GM8BDKO0ahRhFVDAIcOwaWXOqagdesK+wPS0qB3b5g9u1Q7jbP3ZrMgZ0GgbDuJDaNq4mX5aCFU9XQRaQScEwd5jCAaNYLdux1Tf2pq0ImCAlizxlEE4cjMhA8+cN4H7zSGmBzIIz8eWaTOPyt45tJnPF/HSA6OHj1KdnY2hw8frmhRkpqaNWvSvHlzqlXzvrDDS4ayOhzbVNYeOA34RFU/9tB3EPAkkAq8pKqPhGnTD3gCqAbsUNW+nqWv4jRq5Izje/Yc21cAwMaNkJcH7dqF7zh6tDNjCMa/0zgGRTBvw7widbaT2IhEdnY29erVo3Xr1rbXpIJQVXbu3El2djZt2rTx3M/LjGAO0EdEjgdmAAuAa4GoI4qIpALPABcA2cD/RGSqqi4LatMAeBYYpKobROQEz5InAcG7iwspglWrnNdIiiDS8tIYdhpv3b+VA0cP8Kdef+KxCx/z3M9IXg4fPmxKoIIRERo1asT27dtj6ufFRyCqehC4CnhKVa8EzvDQ7yxgtapmqWoe8AZwRUibXwPvqeoGAFXd5l30qk/EMBN+RdC+ffiOkZaXxrDT+JUfXiHfl8/w7sM99zEMUwIVT0n+Bp4UgYj0wpkB+M1BXmYSJwMbg8rZbl0w7YHjRWSWiCyMFN5aREaIyAIRWRCrpqvMRFUEtWvDSSeF7zhunHM+mBh2GqsqL33/En1a9uG0xqfFJrRhGJUOL4rgDuBe4H1VXSoibYEvPfQLp5ZCo6OlAT2AS4GLgPtdP0ThTqovqGpPVe3ZpEkTD7euGkSMQLpqFZx66rFdxKEMHeosJ23V6lib3/7Ws39g9vrZrNq1ymYDRqVk3LhxdOzYkS5dupCens63335bIXIsXryYadOmBcpTp07lkUccN+mNN97IO++8U6TPrFmzuOyyy8pNRj9eYg3NBmYHlbOA2zxcOxtoEVRuDmwO02aHqh4ADojIHJww1z97uH6VJ2LguVWroHPn6J2HDnUOVUch5OR4umdObg7XvH0N9arXI+OMjNiFNgwPNGsGW7cWrW/atHT5ub/++ms++ugjvv/+e2rUqMGOHTvIy8sr+QVLweLFi1mwYAGXXHIJAIMHD2bw4MRM5RJxRiAiH4rI1EiHh2v/D2gnIm1EpDowBAjt918cR3SaiNQGzgaWl/TDVDWOO87ZElBoRpCfD1lZkR3FoYhARgZ89hns21ds89EzR7P94HZa1W9F7Wq1i21vGCUhnBKIVu+VnJwcGjduTA030Xfjxo056aSTaN26NTvcJ6oFCxbQr18/AGbPnk16ejrp6el069aN3NxcAB599FE6d+5M165dueeeewBYs2YNgwYNokePHvTp04cVK1YAztP9yJEj6dOnD+3bt+ejjz4iLy+Pv//977z55pukp6fz5ptvMnnyZP74x2Pbr6ZPn16oTygHDhzgpptu4swzz6Rbt27897//BWDp0qWcddZZpKen06VLF1b5fYalINqMwL9U5CqgGTDFLV8HrCvuwqqa7246+wxn+ehE17Q00j0/QVWXuxFNfwR8OEtMl5Tok1RBRKBhwxBFsG6dowy8KgJwFMHjjzshJ37964jNghPTr9q1yhLTGyXmjjtg8eKS9XXH6CKkp8MTT0Tve+GFF/LQQw/Rvn17Bg4cyLXXXkvfvpFXpD/22GM888wz9O7dm/3791OzZk0++eQTPvjgA7799ltq167Nrl27ABgxYgQTJkygXbt2fPvtt4waNYqZM2cCsG7dOmbPns2aNWs4//zzWb16NQ899BALFizg6aefBmDy5MmF7h2uTzDjxo2jf//+TJw4kT179nDWWWcxcOBAJkyYwO23387QoUPJy8ujoAxiikVUBK5JCBEZq6rnBZ360DXhFIuqTgOmhdRNCCn/E/inZ4mTjCLxhopbOhqOc85xHMvvvBNVEYydMzaQitKfmN42jhmVibp167Jw4ULmzp3Ll19+ybXXXhuwy4ejd+/e3HXXXQwdOpSrrrqK5s2bM336dIYNG0Ztd8FFw4YN2b9/P1999RVXX30sDPuRI0cC76+55hpSUlJo164dbdu2DcwWolFcn88//5ypU6fy2GPOM/nhw4fZsGEDvXr1Yty4cWRnZ3PVVVfRLpaxIAJeVv80EZG2rm8AEWkDJI/HtoIpE0WQkgK/+hW8+CLs3w916xZpkpObw8RFE1HXn+8PJ3F/3/ttVmDETHFP7tFWOM6aVbp7p6am0q9fP/r160fnzp155ZVXSEtLw+fzARTa+XzPPfdw6aWXMm3aNM455xymT5+OqhZZgunz+WjQoAGLI0xzQtt7WcJZXB9V5d1336VDhw6F6k8//XTOPvtsPv74Yy666CJeeukl+vfvX+z9ouFl1dCdwCx3iecsnBVDd5TqroZnwiqCunUdr1osZGTA4cMwbVrY08GzAT8WZM6obKxcubKQzXzx4sW0atWK1q1bs3DhQgDefffdwPk1a9bQuXNn7r77bnr27MmKFSu48MILmThxIgfd5E+7du3iuOOOo02bNrz99tuAM0j/8MMPgeu8/fbb+Hw+1qxZQ1ZWFh06dKBevXoBn0M4wvUJ5qKLLuKpp55C3VS0ixYtAiArK4u2bdty2223MXjwYH78sfQpYopVBKr6KdAOuN09OqjqZ6W+s+GJIhFIV61yNpLFummkd2+oVw+GDQsbiO7r7K8p0MK2RgsnYcSLSM8xsT7fhLJ//35uuOEGzjjjDLp06cKyZct44IEHGDNmDLfffjt9+vQhNShw1xNPPEGnTp3o2rUrtWrV4uKLL2bQoEEMHjyYnj17kp6eHjDNZGZm8vLLL9O1a1c6duwYcN4CdOjQgb59+3LxxRczYcIEatasyfnnn8+yZcsCzuJQwvUJ5v777+fo0aN06dKFTp06cf/99wPw5ptv0qlTJ9LT01mxYgXXXx92+1VsqGrUA6gN/A140S23Ay4rrl+8jh49emgy8de/qlavrurzuRVt26pee23sF5oyRTUtTdVZUOoctWs79aq6aucq5QH0X1/9q+yEN5KKZcuWVbQIFcINN9ygb7/9dkWLUYhwfwtggUYYV72YhiYBeUAvt5wNPFx6FWR4oVEjJ77cgQM4b9ati80/4Gf0aGe1UTD+QHTAm0ucJ5ZrOl5TOoENw6h0eHEWn6Kq14rIdQCqekgsoEi5ERxmou7htU7+gZIogmIC0b2x9A3ObXkuzY9rXkJJDSM5CV0WWhnxMiPIE5FauOEhROQU4Ej0LkZZUSje0M/uhuuSKIIogeiWbFvCkm1LGNJxSIlkNAyjcuNFEYwBPgVaiEgmTijqv8ZVKiNAIUVQkqWjfsIFonPzHL+55E1SJMVCShhGkhLVNCQiKcDxOLuLz8EJJHe7qkZKoGiUMUUUQYMGIckJPOIPOHfffY45qE4deP55Nl9+Po89NZzeLXrTtG4pl2wYhlEpiTojUFUf8EdV3amqH6vqR6YEypciiqBdu9iXjvoZOtRJWzlsmBPEKCOD2z+9ncP5h0mV1OL7G4ZRJfFiGvpCRP4sIi1EpKH/iLtkBuDEGgI3AqlfEZSWjAzYu5ddU9/igxUfAPDNpm8sKb1R6UlNTSU9PZ1OnTpx+eWXs2fPnooWKWYeeOCBwN6F8sKLIrgJ+ANOysqF7rEgnkIZx6hWDerXh71bDzu5iiNlJYuFAQOgfn1WPD+OfHWWlPrUZ7uIjXInJzeHvpP7ltlDSK1atVi8eDFLliyhYcOGPPNMYsTKUtVAiItExMvO4jZhjrblIZzhcGP1TO6f1NbZBvb004V2BJeIGjU4ePEFnD5/JdXcrQX+2EI2KzDKk7FzxjJvw7y4PIT06tWLTZs2AZFDSG/dupUrr7ySrl270rVrV776ytlJ/+9//5tOnTrRqVMnnnADJ9199908++yzges/8MAD/Otf/wLgn//8J2eeeSZdunRhzJgxgBNd9PTTT2fUqFF0796djRs3hm0HTqTRDh06MHDgQFauXBmoHz9+fGCX9JAh8VvVV+w+AjdPwF1AS1UdISLtcMJMFA2gbZQ9mZn8v50jqOVz4p6wYweMGOG895hxLByvtDvALYfh/HXw+alOnT+2kEUcNUrLHZ/eweIti6O2OZJ/hO82f4dPfUxYOIFFWxZRPbV6xPbpzdJ5YtATnu5fUFDAjBkz+N3vfgdEDiF922230bdvX95//30KCgrYv38/CxcuZNKkSXz77beoKmeffTZ9+/ZlyJAh3HHHHYwaNQqAt956i08//ZTPP/+cVatW8d1336GqDB48mDlz5tCyZUtWrlzJpEmTePbZZyO2q1OnDm+88QaLFi0iPz+f7t2706NHDwAeeeQR1q5dS40aNeJq5vKyoWwSjjnoF245G3gbMEVQHowefUwJ+PHvCC6FIph8wiaGVoeMZccUgcUWMsqT9XvXBwKqqSrr96ynXaPS+cAOHTpEeno669ato0ePHlxwwQVRQ0jPnDmT//znP4DjX6hfvz7z5s3jyiuvpE6dOgBcddVVzJ07l9tuu41t27axefNmtm/fzvHHH0/Lli0ZP348n3/+Od26dQOceEerVq2iZcuWtGrVinPOOQdwwkqHa5ebm8uVV14ZCHsdnMWsS5cuDB06lF/+8pf88pe/LNV3Ew3bWZzoFLMjuKR8fNMMfvq/JvxukXDzIpwNZ+PGlUq5GIaf4p7cc3JzaDu+bSDsuaLsPrybNzLeKFXYc7+PYO/evVx22WU888wz3HjjjVFDSIfiV07hyMjI4J133mHLli0BU42qcu+99/L73/++UNt169YFlEm0dk888UTEsNUff/wxc+bMYerUqYwdO5alS5eSluZl2I4N21mc6ETZEVwalj35N3rkQIo/BN369Y7JqbT+B8PwwNg5Y/FpYedpWYY9r1+/PuPHj+exxx6jVq1aEUNIDxgwgOeee865f0EB+/bt47zzzuODDz7g4MGDHDhwgPfff58+ffoAMGTIEN544w3eeecdMjKcDZgXXXQREydOZP/+/QBs2rSJbdu2FZEpUrvzzjuP999/n0OHDpGbm8uHH34IODkQNm7cyPnnn8+jjz7Knj17An3LGi+q5QEK7yzuDdwYF2mMoowbR96wEVQ/GmQecncEl4bT//0qNUMz3JWByckwvPB19tfkFRROKl/Wpslu3brRtWtX3njjDTIzM7nlllt4+OGHOXr0KEOGDKFr1648+eSTjBgxgpdffpnU1FSee+45evXqxY033shZZ50FwPDhwwPmnI4dO5Kbm8vJJ5/MiSeeCDjpMZcvX06vXk5czrp16zJlypRC4a6jtevevTvXXnst6enptGrVKqB0CgoK+M1vfsPevXtRVe68804aNGhQZt9PMBJtGhRoJNKIYzuLv6nITWU9e/bUBQuSa/Xq5zdm0u+VG6lGPtKqValNOHkFeaSl1Qg/HRRxAtsZRowsX76c008/vaLFMAj/txCRharaM1z7Yk1DIjIVuBCYZTuLK4bdFw/lAHXZdd0fnDDUpXxin7N+DhvqRzhZSpOTYRiVDy8+gn8BfYBlIvK2iGSISM3iOhllR+N6RziePeTWKptYQB+u/JAxF6ShtWsVPlEGJifDMCofXjaUzVbVUUBb4AXgGqCoJ8SIG03F+bp31yh9EnlV5cOfP2TXVYOQF16EVq2OnXzySfMPGEYS4mVGgLtq6FfASOBM4JV4CmUUpnG+s9t3e0rpZwTLti9j7Z61XN7+cmfQX7cOvv7aOZlqgecMIxnx4iN4E1gO9AeewdlXcGu8BTOO0eDIVgByfKVXBK/99BoAZ5505rHKs8+GU06xpaOGkaR4zVl8iqqOVNWZbmhqoxypsduZEWTnl940NHHRRABe+v6lY5Ui8JvfwMyZ4MZmMQwjefCiCGYAfxCRd9zjVhGpFm/BjGPINmdGsP5w6WYES7ctZcsBR6kUCTA3dKizsez110t1D8OoSOrWrVtsm7lz59KxY0fS09M5dOhQTNf/4IMPWLZsWaD897//nenTp8csZ6LhRRE8B/QAnnWP7m6dUV5s3cq+lPps3Vu6xVp3fXZX4H2RXZzt2jkmoilTSnUPw/BMZia0bg0pKc5rOZkmMzMz+fOf/8zixYupVatW8R2CCFUEDz30EAMHDixrEcsdL4rgTFW9wTULzVTVYTgOY6O82LKFvTWaOlnKSkhObg4z1s4IlMOGne7QAX74odz/YxpJSGamE9Jk/fq4hDiZNWsW/fr1IyMjg9NOO42hQ4eiqrz00ku89dZbPPTQQwx1V8hFCg39n//8hy5dutC1a1d++9vf8tVXXzF16lT+8pe/kJ6ezpo1a7jxxht55513AJgxYwbdunWjc+fO3HTTTYHAdq1bt2bHDmf71YIFC+jXrx8As2fPJj09nfT0dLp160Zubm6ZfPaS4CXERIGInKKqawBEpC0QGpzAiBPNmsGbW7ciNGX+/GNZKps2hS0xpA4YO2csBVr4z1Yo7HRmJrixWAr9xwRbUmrEzh13QLQgb998A0dCQpYdPAi/+x28+GL4Punp4OYG8MKiRYtYunQpJ510Er1792b+/PkMHz6cefPmcdlll5GRkRExNHSjRo0YN24c8+fPp3HjxuzatYuGDRsyePDgQN9gDh8+zI033siMGTNo3749119/Pc899xx33HFHRPkee+wxnnnmGXr37s3+/fupWbPitmd5mRH8BfhSRGaJyGxgJvCn+Ipl+Nm6FZqylS00K1IfC7PWzSpSVyi2y+jREGov9cceMoyyJlQJFFdfAs466yyaN29OSkpKIDR1KMGhobt3786KFStYtWoVM2fOJCMjg8aNGwPQ0J8zNgIrV66kTZs2tHczCN5www3MmTMnap/evXtz1113MX78ePbs2ROXqKJeKfbOqjrDn4wGJ9bQClW16KPlSDO28AUXlOoad55zJyM+GsHyPyzntManFW0Qp3DXRpJS3JN769bOrDOUVq1g1qwyEaFGjRqB96mpqeTn5xdpEyk09Pjx4yOGhg5HtJhtaWlpgTSVhw8fDtTfc889XHrppUybNo1zzjmH6dOnc9ppYf5vlgNe9hH8Aailqj+q6g9AbREZ5eXiIjJIRFaKyGoRuSdKuzNFpEBEMiK1SVZqcJgG7C0yI4iV6Wunc3K9k+nQqEP4BnEKd20YYRk3zglpEkwFhDiJFBp6wIABvPXWW+x0HXO7du0CoF69emFt+aeddhrr1q1j9erVALz66qv07dsXcHwECxcuBODdd98N9FmzZg2dO3fm7rvvpmfPnoH0mRWBF9PQzaq6x19Q1d3AzcV1EpFUnA1oFwNnANeJyBkR2v0f8JlHmZOKpjg2oK2UfOmoT33MyJrBgLYDIj/lhPuPWauWxR4y4sPQofDCC84MQMR5feGFcvdHXXjhhfz617+mV69edO7cmYyMDHJzc+nYsSOjR4+mb9++dO3albvuclbcDRkyhH/+859069aNNWvWBK5Ts2ZNJk2axNVXX03nzp1JSUlh5MiRAIwZM4bbb7+dPn36FApN/cQTT9CpUye6du1KrVq1uPjii8v1sxdCVaMewI+44ardciqw1EO/XsBnQeV7gXvDtLsD+AMwGcgo7ro9evTQZOJMvlUFvYypbgaZY4dXFuUsUh5A/7P4P9EbTpmi2qqVqohzg1/9qlSyG8nFsmXLKloEwyXc3wJYoBHGVS8zgs+At0RkgIj0B17HSVRTHCcDG4PK2W5dABE5GbgSmBDtQiIyQkQWiMiC7du3e7h11aFDA2dGEGoaahrDBGF6lrPhZUDbAdEb+mMP+XwwcCB89x0U2AIxw6jqeFEEd+OsFLoF58l9BvBXD/3C2SBCPSpPAHeratTRRlVfUNWeqtqzSZMmHm5ddXj1UWeN6Mwlzsj/6KPOfCCWpaPTs6ZzRpMzOKneSd47jRwJGzfCp150vmEYlRkvq4Z8IjIZmKmqK2O4djbQIqjcHNgc0qYn8IZrt24MXCIi+ar6QQz3qdq460TrntKUWrViXzZ6JP8Ic9bP4ebuxbp1CjN4sLOJYcIEuPTS2PoaSYuqxrTaxih71EPWyVC8rBoaDCzGNQeJSLqbtaw4/ge0E5E2IlIdGAIU6qeqbVS1taq2Bt4BRpkSCGHLFmjQAKlZg2bNYpsJgJMb9lD+oeLNQqFUq+Zs7pk2zZkZGEYx1KxZk507d5ZoIDLKBlVl586dMW9O87KDYQxwFjDLvdFiEWntQaB8Efkjjo8hFZioqktFZKR7PqpfwHDZujXgECiJIpieNZ1USaVvq76x3/vmm51VQ507w759zlLSUuZLNqouzZs3Jzs7m2Tz4yUaNWvWpHnz5jH18aII8lV1b0mme6o6DZgWUhdWAajqjTHfIBnYutXRADgvP//svWtObg7jvx1PerN06teMlKQ4CvPmOXGH9u51yhZ2wohCtWrVaNOmTUWLYZQAL87iJSLyayBVRNqJyFPAV3GWy/CzZUuJZwR/m/k3cvNKEchq9GhnBVEwFnbCMKocXhTBrUBH4AjO0tG9wO3xFMoIImRGsHMn5OUV3y0nN4cpPzkhpZdsW1I4yqhXLOyEYSQFXpLXH1TV0ap6pqr2BKYAT8dfNINDhxzbfNCMAGDbtuK7jp0zlnyfE1tF0cK5B7xiYScMIymIqAhEpIuIfC4iS0RkrIg0FZF3genAskj9jDLEv1Y0RBEUZx7Kyc1h0uJJ+NysomFzD3ghQeLBGIYRX6LNCF4EXgN+BewAvgeygFNV9fFykM3wK4Ig0xAUrwjGzhkbUAJ+imQk80JwPBg/Dz5ojmLDqGJEUwQ1VHWyqq5U1ScBH3CPqh6O0scoS/wjfowzgq+zvyavoLAjoVDugVjwh53YuhVq1oxt2ZJhGJWCaMtHa4pIN46FitgPdBF3Hamqfh9v4ZKekBmBP75QcYpg0e8Xcd+M+/jnV/9kz917qFO9TullOeEEGDYMXn4ZHnromFYyDKPSE21GkAP8G/iXe2wJKj8Wf9GMwIh/wgkA1KgBxx/vbQnp3A1z6XFij7JRAn7uustZstShg+U1NowqRMQZgaqeX56CGGHYutUZ+atXD1R52UtwOP8w3236jtvOuq1s5fn2W0hNdVYygW0wM4wqgpd9BEZFEbSHwI+XpPXfbfqOvII8zmt1XtnKM3p00bDUtsHMMCo9pggSmaBdxX68zAjmrHeSZvdu2bts5bENZoZRJTFFkMiEmRF4UQRzN8yl8wmdaVirYdnKYxvMDKNKErMiEJETRaRGPIQxQogwIzhwANxc20XI9+Xz1cav6NOyT9nLY3mNDaNKUpIZwavAChGxlUPx5OBBZ7QPowggcoKaxVsWsz9vf9n7B6BownGAiy4yR7FhVHJiVgSqOhBoC0wqe3GMACF7CPwUt6ls7vq5APRpFYcZARTOa3zNNfDFF7GnTTMMI6HwkqHsFL8pSET6ichtQH1VXRp36ZKZkF3FfopTBHM2zOGU40+JLT9xSRk7Fg4fhn/8I/73MgwjbniZEbwLFIjIqcDLQBucGERGPCnBjEBVmbdhXvxmA6G0bw99+sD48bbBzDAqMV4UgU9V84ErgSdU9U7gxPiKZYRGHvXTuLEz5oZTBHPXz2XHwR10OaFLOQiIM+h/+63zXvXYBjNTBoZRqfCiCI6KyHXADcBHbl21+IlkAEXCS/hJTXWqwimC+7+8H4CFOQvjLZ3D6NFOzoRgbIOZYVQ6vCiCYUAvYJyqrhWRNjjJaYx4kZkJj7mLstq1K/KEHW4vQU5uDvM3zgfgveXvlSwjWazYBjPDqBJ4yVC2DLgbJx8BqrpWVR+Jt2BJS2amY17xbxQIY24JpwjGzhlLgTrhH0qUe6Ak2AYzw6gSeFk1dDmwGPjULaeLyNQ4y5W8jB7tmFeCCTG3hCqCnNwcJi6aGCiXOCNZrITbYAZw333xva9hGGWKF9PQA8BZwB4AVV2Ms3LIiAcezC3Nmjm+ZJ+bhCx4NuCnXGYFoRvMTnTXECxeHN/7GoZRpnhRBPmqujekTuMhjIEnc0vTpnD0KOze7ZS/zv46kKjeT4kzksVK8AazzZvhttvgueccpWBLSg2jUhAtQ5mfJSLyayBVRNoBtwHlMMIkKePGwfDhzkYtPyEJ44P3EjRq5GQk6/58dxrWasj066eXs8AhdO58TDiwnAWGUQnwMiO4FegIHAFeB/YBd8RRpuRm6FC4/XbnvYhjdnnhhUKDaOimstwjufyw9Qd6tyjjsNMl4eGHi9bZklLDSGiKnRGo6kFgtHsY5UHr1s7r+vXQokWR06GK4Jvsb/Cpj3Nbnls+8kXDlpQaRqUjoiIQkQ+J4gtQ1cFxkciArCwnQfHJJ4c9HRqBdP7G+aRICuc0P6ecBIxCy5aOAgtXbxhGQhJtRmBhpiuKNWugTRvH2RqG+vUdPeGfEczbMI8uTbtQr0a9chQyAuPGOT6B4CWw1apZzgLDSGAi+ghUdbb/AL4GdgO7gK/dOiNeZGVB27YRT4sc20uQ78vnm+xvOLdFApiFoOiS0jp1nCVOf/6zrSIyjATFy4ayS4E1wHjgaWC1iFzs5eIiMkhEVorIahG5J8z5oSLyo3t8JSJdY/0AVQ5VZ0ZwyilRm/kVwQ9bfuDA0QNln5+4NAQvKX38cUchbNligekMI0HxsmroX8D5qtpPVfsC5wOPF9dJRFKBZ4CLgTOA60TkjJBma4G+qtoFGAu8EIvwVZKdOyE3N+qMAI4pAn98oYRYMRSOceMcBRCMrSIyjITCiyLYpqqrg8pZwDYP/c4CVqtqlqrmAW8AVwQ3UNWvVNXdFsU3QHMP163arFnjvHqcEczfOJ+W9VvSon7R1UUJga0iMoyEJ9qqoavct0tFZBrwFs4qoquB/3m49snAxqByNnB2lPa/Az6JIMsIYARAy6q++iQry3n1MCPYvkOZt34efVv3LQfBSoitIjKMhCfajOBy96gJbAX6Av2A7cDxHq4tYerCLkcVkfNxFMHd4c6r6guq2lNVezZp0sTDrSsx/hlBm8jhnJo1gwcfBOqvZ/P+zbz+f70DDuSEI1Jgun37zHlsGAlCxBmBqg4r5bWzgWB7RXNgc2gjEekCvARcrKo7S3nPyk9WlhOnJ9zg6RLIFd9ynvO64dzC9YmEf0f06NGOOei442Dv3mOBkiwEhWFUOBFnBCLyV/f1KREZH3p4uPb/gHYi0kZEqgNDgELhq0WkJfAe8FtV/bnkH6MK4WHFUIC2n4MvFQ40jq9MpSV4FVGDBkXPm/PYMCqUaBvKlruvC0pyYVXNF5E/Ap8BqcBEVV0qIiPd8xOAvwONgGdFBJxIpz1Lcr8qQ1YW9O/vrW37j0EK4Lx/wLRn4itXWWHOY8NIOKKZhj503x5U1beDz4nI1V4urqrTgGkhdROC3g8HhnuWtqpz+DBs2uRtRtBoBdTa5Xhiuk2COffD/kR0EoRgzmPDSDi8LB+912OdUVrWrXPW3BezYgiAQXcec8dLAZxXDqkpy4JIzuPTT3ccx+ZANoxyJ9ry0YuBS4CTQ3wCxwH54XsZpcLjHoLGbXLYccoXxyrS8qDbJBovux9I8FlBqPO4eXM4cgQ+/fRYG3MgG0a5Em1GsBnHP3AYWBh0TAUuir9oSYhfERQzI7j6qbGQUjg1ZfWaBVzzdCWZFQQ7jzdsgOrVi7YxB7JhlBvRfAQ/AD+IyGs4RojTcPYBrHR3ChtlTVaWE6TthBOiNpuzfk6RunJLTRkPNm0KX28OZMMoF7z4CC6ghEHnjBhZs8aZDUi4vXjHuOdcJ37fgpsXcsKzyg1rFR2jLPr9ovKQsuyJ5ChOSTGfgWGUA14Uwb8pQdA5owRkZXlaMTRj7Qwa1mpItxPT6dEDFi4sB9niSSQHckGBRSw1jHIgnkHnjFjw+TwpAlVletZ0+rfpT4qk0KMHLFtWOA9MpSM0h0FqatE25jMwjLjhRREsFZFpInKjiNwAfAj8T0SuCgpMZ5SWLVucfQTFOIpX7VpF9r5sBrQZAECPHo4O+fHH8hAyjgQ7kH2+8G3MZ2AYccGLIggXdK4hTkC6y+ImWbLhceno9KzpAAxsOxBwFAFUAfNQMOYzMIxyJVqICaBMgs8ZXvAYfnrG2hm0rN+SU453FEbz5tC4cRVTBOHyHoPjMwDbZ2AYZYyXVJXtRWSGiCxxy11E5G/xFy3JWLPGedpt1SpikwJfAV+u/ZKBbQbixmZChKrhMA7GfAaGUa54MQ29iBNS4iiAqv6IE0nUKEuysqBFi/Cbq1wWbVnE7sO7GdB2QKH6Hj1g6VI4dCjeQpYjXnwG69dbWArDKAO8KILaqvpdSJ2FmChrPISfnpE1AyDgKPbTo4djNfnpp7hJV7FEC0i3fr0tMTWMUuJFEewQkVNws4uJSAaQE1epko3MTPjuO5g5M+qT7fS10+l0Qiea1m1aqL5KOoyDibTPIBQzFxlGifCiCP4APA+cJiKbgDuAW+IpVFKRmQk333zM/BHhyXbt7rXMyJpBr5N7FblEy5bQqFEVVgShPoMofhTWrzdTkWHESLGKQFWzVHUg0AQ4TVXPVdV1cZcsWRg9uqhxP8yT7e2f3o6ibMotGpenSjqMQwn2GaxbF10ZmKnIMGLCy6qhf4hIA1U9oKq5InK8iDxcHsIlBR4yduXk5jBtlZPfZ+a6mWzZv6VI8+7dYckSZ09aUuDFXHTwINx+uzmUDaMYvJiGLlbVPf6Cqu7GyVNglAUtWoSvD3KQjpk1hgJ11tD71MfY2YXDTTdrBo88Avn5UKuWM0MQceqrLKHmokjs3GkOZcMoBi+KIFVEavgLIlILqBGlvRELN99ctK52beeJF2c2MHnx5MCpvII8Ji2eVGhWsHVr+EtHqq8yBJuLopmKgjl4EG64wWYIhhGEF0UwBZghIr8TkZuAL4BX4itWErF3r/NE27z5MUfoCy8EdsyOnTOWfF/h1boFWlBkVpD0eF1ZBBbV1DBC8OIsfhR4GDgd6AiMdeuM0uLzwVtvwSWXwMaNxxyhQWETZq6diTordwNU6iQ08SLcyqJGjYrvZzMEw/DkLK4DfK6qfwZeAGqISLW4S5YMfPON4xS+9tqITc5rdR610mqx++7d6BgNHJU2CU08CV1Z9OST3mYJoTOEUaPMwWwkFV5MQ3OAmiJyMjAdGAZMjqdQScMbb0CNGnDFFWFP7z28l8yfMrmu03U0qNmgfGWrCniJWRTKwYMwYUJhB/OwYU5kP1MMRhXFiyIQVT0IXAU8papXAmfEV6wkoKAA3n4bLr0UjjsubJNXf3yVg0cPcsuZ0ffvNW0aW31SETxLeOUVbzMELWyK4+hRZ/WR+RWMKoonRSAivYChwMduXbHhq41imDPHSUYTwSy0ed9m7p5+N12bdqXnST2jXmrLFmeM8h+33uo8/M6dGw/BKzElmSGEw/wKRhXDiyK4Ayf66PuqulRE2gJfxlWqqkxmpjN49O/vDEb794dtNmraKA4ePUjDWg1jvsUbbzgTjvbtj+0pqPL7CrxS3Awh2p6EYIrzK5ifwahEiIZOgxOcnj176oIFCypajJKRmVk04Urt2oWWi4Kzd6D5483xqY9aabXIuj2LZnW9j+LRxrJK9ueOP5mZTjiPDRucTXyXXOIoiLJOAl2tmmMC3LXLuc+4cZZUxyhXRGShqoY1L3hZNfSliMwMPcpezCRg9OiiA0yYuEJXvnklPnWC0NmegTgTutLo2WcLm48aNYqaI8Iz4fwMNoswEoRiZwQi0iOoWBP4FZCvqn+Np2CRqNQzgpSU8I/kIoHoo/9v7v/jvpn3FTod66zAq3UDHIfylqKhi4xggmcNKSnHUmbGm3CzCCg8g7GZheGRUs0IVHVh0DFfVe8Czi5zKasqfp9ASkrkEbplS3Jyczjt6dO4b+Z9CIXbxXNWUOXDUJQFZeVXiJXQWcSwYXDTTdGXtoabVQT/Bm3mYYTBi2moYdDRWEQuAsztGIng/3SNGxf+jxsu5aIbV+im/97Eyp0rSUtJK/edxMEO5dDDHMwhhNvBPHKk9/AWpeHoUcjLK1oXrCyee66ooghVHuHahO6TKE55lFThlOS6RvxR1agHsBbIcl9XAZ8D5xbXz+07CFgJrAbuCXNegPHu+R+B7sVds0ePHhozU6bo0RYnawHo0RbNVadMKVp3yy3Ryx77aO3awSs5wx75KaIFoDua1NXvHr1DL8u8THkA5QG05sM1NSc3J/bPGETTpsWKYEcZH9cxRdfSSgsQXUsrfYpbAuVtNNLDVC/UoQCpeKEjHIepVkReX0ib0PIhD328tAkt76d2oe8y9LsNV76OKVH/HlWhz3VM0aZNYxsXgAWqEcbqSCdKewCpwBqgLVAd+AE4I6TNJcAnrkI4B/i2uOvGrAimTCk6OFerplq98A+y2MNLH/H2n7sAAgM/D6A1x9bUlAdTlAfQ6mOr66iPRsX2GT2QAONLUh/h/vPvp/DvMnQQtMM5CkLKxSmPeCmyROqzn9p6HVNiHAMiK4KIpiEROVNEmgWVrxeR/4rIeBHxsrj9LGC1OhnO8oA3gCtC2lwB/MeV8xuggYic6OHa3gm3UifcFLs4vPRR9XSpDfWd1xRSGNhmIAiBVULhwkwblZ/XGUob1pGKjzas41ae5WZeYB2t8CGsoxXPcEugvJ1GHKHwaqXDVCtSlwyEDlKhHpnQcg2OUoO8qG1CyzUrWZ86HOQflF1+7mg+gufBubuInAc8AvwH2IsTfK44TgY2BpWz3bpY2yAiI0RkgYgs2L59u4dbBxEpA1gFcaAa3DfAee/Dx5frvgwoAT/xcA5buInEI5xy8JdPYAfDmFhIUdzEpEJ14ZRF6KNIOOXh7XHFSHRaUnZjWzRFkKqqu9z31wIvqOq7qno/cKqHa4dbShH6G/TSBlV9QVV7qmrPJk2aeLj1MfKbnxRT+9IS6g4+nALbazn16+rDzZfD612OnS/QAvIKCmv7eDiHQ8NQGIlPqKJ4naGF6sIpi+BZRTjlUdKZR+hPpiQKpyTX9YUdIgyADbQsvpFXItmMgCVAmvt+BXBe8LlI/YLa9AI+CyrfC9wb0uZ54Lqg8krgxGjXjdVH8PKfBuj+aiE2txT0cGpsdkcvffZXQ5/uia6t79g119ZHr7uqsD8g3JE+IT2mz1QWmEPZDv8RzhEZDydoSa7rxZeSSLb7yuojCFvp9GE0MB/4L7CIY5vPTgXmR+oX1D8NZ7VRG445izuGtLmUws7i74q7bqyKIH1Cul53VdHBObTuqZ7Ry176eBn0K2rgLwmmLOxIhKMyreaprKuGou4sFpFzgBNxEtMccOvaA3VV9fviZhsicgnwBM4KoomqOk5ERrozkQkiIsDTOMtMDwLDVDXqtuFKvbPYMAyjgoi2s9iCzhmGYSQBpQoxYRiGYVRtTBEYhmEkOaYIDMMwkhxTBIZhGElOpXMWi8h2YH0JuzcGdpShOPGmMslbmWSFyiVvZZIVKpe8lUlWKJ28rVQ17I7cSqcISoOILIjkNU9EKpO8lUlWqFzyViZZoXLJW5lkhfjJa6YhwzCMJMcUgWEYRpKTbIrAS9TURKIyyVuZZIXKJW9lkhUql7yVSVaIk7xJ5SMwDMMwipJsMwLDMAwjBFMEhmEYSU7SKAIRGSQiK0VktYjcU9HyhCIiE0Vkm4gsCaprKCJfiMgq9/X4ipTRj4i0EJEvRWS5iCwVkdvd+oSTV0Rqish3IvKDK+uDiSqrHxFJFZFFIvKRW05kWdeJyE8islhEFrh1iSxvAxF5R0RWuL/fXokor4h0cL9T/7FPRO6Il6xJoQhEJBV4BrgYOAO4TkTOqFipijAZJxx3MPcAM1S1HTDDLScC+cCfVPV0nDwSf3C/z0SU9wjQX1W7AunAIDe8eiLK6ud2YHlQOZFlBThfVdOD1rcnsrxPAp+q6mlAV5zvOeHkVdWV7neaDvTACdP/PvGSNVKigqp04CFbWiIcQGuCsr8RlLENJy/EyoqWMYLc/wUuSHR5gdrA98DZiSor0Nz9D94f+CjRfwfAOqBxSF1CygscB6zFXSST6PIGyXchbjKweMmaFDMC4GRgY1A5261LdJqqag6A+3pCBctTBBFpDXQDviVB5XVNLYuBbcAXqpqwsuIkcvorhdNfJ6qsAAp8LiILRWSEW5eo8rYFtgOTXNPbSyJSh8SV188Q4HX3fVxkTRZFEC4Dtq2bLSUiUhd4F7hDVfdVtDyRUNUCdabYzYGzRKRTBYsUFhG5DNimqgsrWpYY6K2q3XHMrn8QkfMqWqAopAHdgedUtRtwgAQwA0VDRKoDg4G343mfZFEE2UCLoHJzYHMFyRILW0XkRAD3dVsFyxNARKrhKIFMVX3PrU5YeQFUdQ8wC8cXk4iy9gYGi8g64A2gv4hMITFlBUBVN7uv23Bs2GeRuPJmA9nujBDgHRzFkKjygqNgv1fVrW45LrImiyL4H9BORNq4GnYIMLWCZfLCVOAG9/0NOLb4CsfNNf0ysFxV/x10KuHkFZEmItLAfV8LGAisIAFlVdV7VbW5qrbG+Y3OVNXfkICyAohIHRGp53+PY8teQoLKq6pbgI0i0sGtGgAsI0HldbmOY2YhiJesFe0IKUeHyyXAz8AaYHRFyxNGvteBHOAozpPL74BGOI7DVe5rw4qW05X1XBzT2o/AYve4JBHlBboAi1xZlwB/d+sTTtYQuftxzFmckLLi2Nx/cI+l/v9XiSqvK1s6sMD9PXwAHJ+o8uIsbtgJ1A+qi4usFmLCMAwjyUkW05BhGIYRAVMEhmEYSY4pAsMwjCTHFIFhGEaSY4rAMAwjyTFFYBgREJFGQdEft4jIJvf9fhF5tqLlM4yywpaPGoYHROQBYL+qPlbRshhGWWMzAsOIERHpF5Qr4AEReUVEPndj818lIo+6Mfo/dUNxICI9RGS2G5ztM3+YAMNIBEwRGEbpOQW4FLgCmAJ8qaqdgUPApa4yeArIUNUewERgXEUJaxihpFW0AIZRBfhEVY+KyE9AKvCpW/8TTo6JDkAn4AsnTBOpOOFEDCMhMEVgGKXnCICq+kTkqB5zvPlw/o8JsFRVe1WUgIYRDTMNGUb8WQk0EZFe4ITwFpGOFSyTYQQwRWAYcUZV84AM4P9E5AecaK2/qFChDCMIWz5qGIaR5NiMwDAMI8kxRWAYhpHkmCIwDMNIckwRGIZhJDmmCAzDMJIcUwSGYRhJjikCwzCMJOf/A5lJqdRWh9IZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# -*- coding: utf-8 -*-\n",
    "\n",
    "###################################\n",
    "### Written by Ilias Soumpasis    #\n",
    "### ilias.soumpasis@ucd.ie (work) #\n",
    "### ilias.soumpasis@gmail.com\t  #\n",
    "###################################\n",
    "\n",
    "import scipy.integrate as spi\n",
    "import numpy as np\n",
    "import pylab as pl\n",
    "\n",
    "beta=1.4247\n",
    "gamma=0.14286\n",
    "TS=1.0\n",
    "ND=70.0\n",
    "S0=1-1e-6\n",
    "I0=1e-6\n",
    "INPUT = (S0, I0, 0.0)\n",
    "\n",
    "def diff_eqs(INP,t):\n",
    "\t'''The main set of equations'''\n",
    "\tY=np.zeros((3))\n",
    "\tV = INP\n",
    "\tY[0] = - beta * V[0] * V[1]\n",
    "\tY[1] = beta * V[0] * V[1] - gamma * V[1]\n",
    "\tY[2] = gamma * V[1]\n",
    "\treturn Y   # For odeint\n",
    "\n",
    "t_start = 0.0; t_end = ND; t_inc = TS\n",
    "t_range = np.arange(t_start, t_end+t_inc, t_inc)\n",
    "RES1 = spi.odeint(diff_eqs,INPUT,t_range)\n",
    "\n",
    "print(RES[:3])\n",
    "\n",
    "#Ploting\n",
    "pl.plot(RES1[:,0], '-bs', label='Susceptibles')  # I change -g to g--  # RES[:,0], '-g',\n",
    "pl.plot(RES1[:,2], '-g^', label='Recovereds')  # RES[:,2], '-k',\n",
    "pl.plot(RES1[:,1], '-ro', label='Infectious')\n",
    "pl.legend(loc=0)\n",
    "pl.title('SIR epidemic without births or deaths')\n",
    "pl.xlabel('Time')\n",
    "pl.ylabel('Susceptibles, Recovereds, and Infectious')\n",
    "#pl.savefig('2.1-SIR-high.png', dpi=900) # This does, too\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-23T12:44:07.431514Z",
     "start_time": "2021-09-23T12:44:07.365317Z"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df = pd.DataFrame(RES, columns = ['Susceptibles', 'Recovereds', 'Infectious'])\n",
    "df.to_excel('Fig2.4-SIR.xlsx', index = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-23T12:56:49.345531Z",
     "start_time": "2021-09-23T12:56:48.199208Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAzEElEQVR4nO3deXjU9bn38fcHCCACVQiyiBD0CG5AUERxqbiDVXyOB+tC69Ze1G7q6aY9PFZrT/p41J7L2tpSatGqqbu11NJK1YJbXUBxAaXKHtkCiLIKyP388f0NTIbJ5DfJTGaSuV/XlWvmt86dIcw9311mhnPOudLVptABOOecKyxPBM45V+I8ETjnXInzROCccyXOE4FzzpU4TwTOOVfiPBG4jCT9l6S7MhxfLOm0PLzueEnTc33fpPv/VdKlGY7fI+m/8/X6mUi6UdL9GY7PlTQqy3vOkPTVpsaWaw39rsV679bGE0ERkXSCpJckfSxpnaQXJR0dHbtM0gtxzs0lM/upmTX7B4iZVZvZGXm8/xgz+z3s+d7mk6QKSSapXWPvYWaHm9mMDK9Rch+AkkZJqil0HC1Vo/8YXW5J6go8CXwdeBhoD5wIfNqUc11paUqCyTVJbc3ss0LH4RrmJYLiMRDAzB4ws8/MbIuZTTezt5p4LpLaSLpO0gJJayU9LKlbdCzxDXWCpOWSVkj6btK1db5dSvqypCXRfSY24nUul7RM0keSrpR0tKS3JK2X9Muke6WWgA6X9Peo9LNK0n+l+T0HRPdpE23fJWl10vH7JV0TPZ8h6auSDgUmASMlbZS0PumW+0r6i6QNkl6RdFDSvY6T9FpUIntN0nFJx+pUl6W8h89Fj+uj1xuZ7t8M6Cjpoei1X5c0NN39o3s/Gv1unwBXAv8FXBDd/82ke/aPSo4bJE2XVB7do2N0/dro/XtNUs90QUk6NHrv1itUUY1NOnaPpF9LmiZpE3BymusHSJoZxfB3oDzl+LEKJd31kt5UUhVY9LfzbnTtQklfi/bvDfwV6BP9zhsl9Ykuay/p3uiauZKGJ93vWkkfRsfmSzq1nn+L1s/M/KcIfoCuwFrg98AYYN+U45cBL8Q5N829rwFeBvoCHYDfAA9ExyoAAx4A9gYGA7XAadHxG4H7o+eHARuBz0f3+V9gR9K5cV5nEtAROAPYCjwB7AfsD6wGTkrz+3YBVgDfja7tAhxTz++6FDgqej4fWAgcmnRsWPR8BvDV1NdKus89wDpgBKHkXA08GB3rBnwEfDk6dlG03T06vjjxnqR5DxPvQ7sM/143AtuBcUAZ8D1gEVCWev+kc/8P4YvdXsmvl3TPGcACwpeIvaLtm6NjXwP+DHQC2gJHAV3TxFUGfEBINO2BU4ANwKCk9+xj4Pgolo5p7vFPwt9NB8Lf0Yak92Z/wt/1WdH1p0fbPaLjXwAOAgScBGwGjoyOjQJq0ryPW6P7tQX+H/BydGwQsAzok/TvclChPwcK9eMlgiJhZp8AJxA+JH4L1Eqamu6bWTbnRr4GTDSzGjP7lPAfZJzqViP82Mw2mdnbwN2ED7dU44Anzey56D7XAzuzfJ2fmNlWM5sObCIkitVm9iHwPDAszeueDaw0s59F124ws1fq+V1nAidJ6hVtPxptDyAk0DfruS6dx83sVTPbQUgEldH+LwDvm9l9ZrbDzB4A3gPOyeLeDZltZo+a2XbCB2dH4Nh6zv2nmT1hZjvNbEuGe95tZv+KznmY3b/PdqA78G8WSpizo7+xVMcCnQkJZJuZPUuookz+W/mTmb0YxbI1+WJJ/YCjgevN7FMze46QgBK+BEwzs2nR9X8HZhE+yDGzv5jZAgtmAtMJVaKZvBDd7zPgPiBRsvqMkIwOk1RmZovNbEED92q1PBEUETN718wuM7O+wBFAH+D2pp4L9Af+GBW31wPvEv4jJCeOZUnPl0T3S9Un+Twz20T4xpbN66xKer4lzXbnNK97AOHbbBwzCd8OP0+ohplB+PZ4EvC8me2s98o9rUx6vjkptj6E9yjZEsI32lxJfp93AjWk/zepc24D6vt97gOeAh5UqB68RVJZmuv7AMtS3sPU3ztTLH2Aj6K/m+TrE/oD5yf+fqK/oROA3gCSxkh6OaoeXE9IEHWqltJI/Z07SmpnZh8QSrA3AqslPZhUnVRyPBEUKTN7j1DUPiIH5y4DxpjZPkk/HaNv4QkHJD3vByxPc58VyedJ6kT4JpnN6zTGMkKVQBwzCd8SR0XPXyBUVZwUbaeT7RS8ywkfWsn6AYnfcxOhmiWhV9LzuK+V/D63IVS3pfs3SXfPrH4fM9tuZj82s8OA4wglsEvSnLocOCDRBhNJ/r0beu0VhHaXvVOuT1gG3Jfy97O3md0sqQPwGHAb0NPM9gGmEaqJGnrdtMzsD2Z2AuHf0oD/yfYerYUngiIh6RBJ35XUN9o+gFDkfrkp50YmAVWS+kfn95B0bso510vqJOlw4HLgoTT3eRQ4W6HranvgJur+DcV5ncZ4Eugl6RpJHSR1kXRMuhPN7H1CyeJLwHNRFccq4D+oPxGsAvpGv1Mc04CBki6W1E7SBYT2kyej43OACyWVRY2T45KurSVUpx3YwGscJem8qFrtGkKPsPr+fVOtAipSPrDrJelkSYMltQU+IVQVpevt8wohyf0g+t1GEarDHozzOma2hFDV82NJ7SWdQN3qtPuBcySdKalt1Ig9Kvo7b0+oyqkFdkgaQ2hnSv6du0v6XMzfeZCkU6IEs5XwN1OyPZw8ERSPDcAxwCtRj4uXgXcIDaRNORfg58BUYLqkDdH5qR+kMwkNgc8At0V1+HWY2Vzgm8AfCN/uPiJUWWTzOlkzsw2EhsNzCEX990nTIyXld1lrZkuTtgW8Uc/5zwJzgZWS1sSIZy3hW/N3CVVjPwDONrPEtdcTSjAfAT8mvF+JazcDVcCLUfVHffX+fwIuYHej9HlRe0Ecj0SPayW9HuP8XoQk/wmhOm8m4UO5DjPbBowldFBYA/wKuCQqkcZ1MeFvYh1wA3Bv0v2XAecSGqNrCSWE7wNtor+BqwhtGx9F95madO17hA4PC6P3taFqng7AzdHvsZLQYWGPnmilQma+ME0pk1TB7h4pOwocjnOuALxE4JxzJc4TgXPOlTivGnLOuRLnJQLnnCtxRTNBVVzl5eVWUVFR6DCcc65FmT179hoz65HuWItLBBUVFcyaNavQYTjnXIsiKXU0/C5eNeSccyXOE4FzzpU4TwTOOVfiPBE451yJ80TgnHMlLm+JQNIUSaslvVPPcUm6Q9IHCksVHpmPOHr1Aqn5fnr1ajgm55wrJvksEdwDjM5wfAxwcPQzAfh1PoJYtarhc3L9ep4cnHMtSd4SQbQM3boMp5wL3BstO/cysI+k3vmKp5BSk4MnBudcMSlkG8H+1F3WroZ6lvqTNEHSLEmzamtrmyW4fGruUopzzmVSyESgNPvSzoBnZpPNbLiZDe/RI+0Iaeecc41UyERQQ911cjOtyeqccy5PCpkIpgKXRL2HjgU+NrMVuX6Rnj1zfcfc8DYD51yxyNukc5IeAEYB5ZJqCOuTlgGY2STCAuBnEdbJ3UxYMD3nVq7Mx13T69WrcfX/3mbgnCukvCUCM7uogeNGWAi91UhNOo1NDM4515xa3DTULUlqYlC65nHnnCswn2LCOedKnCcC55wrcZ4ImlF9PZiKtWeTc640eCJoRitXgln42bgROnWCK69s3p5NzjmXyhNBgey9N5xzDjz2GOzYUehonHOlzBNBAX3xi1BbCzNmFDoS51wp80RQQGPGhJLBww8XOhLnXCnzRFBAe+0FY8fC44/D9u2FjsY5V6o8ERTYBRfA2rXwj38UOhLnXKnyRFBgX/taeDzzTJ+IzjlXGJ4ICqy+uYiaZY6i6mqoqIA2bcLjN75Rd7u6uuFzivmaYorFrymuWFr6NdXV6f5HN5rC3G8tx/Dhw23WrFmFDiNnMs0/lPN/mupqmDgRli6Fbt1gwwbYtq3+88vKQoCZzinma4opFr+muGJp6dd06gSTJ8P48bFvI2m2mQ1Pe8wTQWE1WyKoroYJE2Dz5hze1DlXMP37w+LFsU/PlAi8aqi1Si1KfutbngSca02WLs3ZrXwa6tYo9dv/kiWFjcc5l3v9+uXsVl4iKLC8TEQ3cWJuvv2XlUH79i33mmKKxa8prlha+jWdOkFVVXb3ycATQYElJqL7+tdD+21iUromTUQXt8hYVgbdu4eGiv79QxD9++/evvtumDKl7r7Uc4r5mmKKxa8prlha+jVZNhQ3xBuLi8Qtt8C118LHH0PXrk282QEHQE3Nnvu7d4fOnUOi6NcvfKPI4R+Tc654eWNxCzBgQHhctKiJNzKD/fbbc3+nTvDzn4deBjt3hkdPAs45PBEUjSYnguReQq+/Dscdl9eipHOu9fBeQ0WiSYkg3RiBOXP8w985F4uXCIpEt27QpUsjE0G6XkKbN4f9zjnXAE8ERUIKpYJGJYL6egnlcMCJc6718kRQRAYMyGrE+G71DSzJ4YAT51zr5YmgiCRKBFn36K2qCmMCkuV4wIlzrvVqMBFIOkhSh+j5KElXSdon75GVoIoK2LQJ1qzJ8sLx42HIEGjXznsJOeeyFqdE8BjwmaR/A34HDAD+kNeoSlSjew6ZwbJlcPHFPkbAOZe1OIlgp5ntAP4duN3M/hPond+wSlOjE8GCBbB6dRg74JxzWYqTCLZLugi4FHgy2leW4XzXSI1OBC+9FB6PPz6n8TjnSkOcRHA5MBKoMrNFkgYA9+c3rNLUuTOUlzcyEXTtCocdlpe4nHOtW4OJwMzmmdlVZvZAtL3IzG6Oc3NJoyXNl/SBpOvSHP+cpD9LelPSXEmXZ/8rtC6NGkvw4oswcmSYXsI557LU4BQTkhYBe3RoNLMDG7iuLXAncDpQA7wmaaqZzUs67ZvAPDM7R1IPYL6kajPLYkHP1mXAAHjjjSwuWL8e5s6FL34xXyE551q5OHMNJU9b2hE4H+gW47oRwAdmthBA0oPAuUByIjCgiyQBnYF1wI4Y9261BgyAJ54InX9ifcF/5ZXQa8gbip1zjRSnamht0s+HZnY7cEqMe+8PLEvaron2JfslcCiwHHgbuNrMdqbeSNIESbMkzaqtrY3x0i1XRQVs2wbLl8e84MUXQ8YYMSKfYTnnWrE4VUNHJm22IZQQusS4t9LsS61iOhOYQ0gsBwF/l/S8mX1S5yKzycBkCAvTxHjtFiu551DfvjEueOklGDo0zFjnnHONEKdq6GdJz3cAi4E4FdI1wAFJ230J3/yTXQ7cbGGZtA+i9ohDgFdj3L9VSk4EJ57YwMk7doSqoUsvzXtczrnWq8FEYGYnN/LerwEHR91NPwQuBC5OOWcpcCrwvKSewCBgYSNfr1VIrCUTq+fQ22/Dxo3ePuCca5I4VUOfA24APh/tmgncZGYfZ7rOzHZI+hbwFNAWmGJmcyVdGR2fBPwEuEfS24SqpGvNLNuZdlqVDh2gT5+YicAHkjnnciBO1dAU4B12Vwd9GbgbOK+hC81sGjAtZd+kpOfLgTPiBlsqYo0lqK6GH/wgPP/85+GnP/X5hZxzjRInERxkZv+RtP1jSXPyFI8jJIIZMzKckLo05dKlYRs8GTjnshanp/oWSSckNiQdD2zJX0huwAD48MPQjTQtX5rSOZdDcUoEXwd+H7UViDDo67J8BlXqKirCgLJly+Cgg9Kc4EtTOudyKE6voTnAUEldo+1PMl/hmiq5C2naRNCvHyxZkn6/c85lqd5EIOlLZna/pO+k7AfAzP43z7GVpF69YNWq8Pz003fv79kTVq6MNqqq4PLLYfv23Sf40pTOuUbK1Eawd/TYJc1P5zzHVbISSSDj/vHj4cwzw3NfmtI510T1lgjM7DfR06fN7MXkY1GDsSukbt3qryJyzrksxOk19IuY+1xzWr48jDxzzrkmytRGMBI4DuiR0k7QlTBS2BXS8uUwaFCho3DOtQKZSgTtCW0B7ajbPvAJMC7/obmMVqzwEoFzLicytRHMBGZKusfMvCK6mfTsmb7BuGfPpI0tW+CjjzwROOdyIk4bwV2S9klsSNpX0lP5C6m0rVwZFhy75ZawvWFD2N7VdRRCaQCgd+9mj8851/rESQTlZrY+sWFmHwH75S0iB0B5eXhMuyBbIhF4icA5lwNxEsFOSbuGrErqT5rF7F1u9egRHtekm5Q7sY6lJwLnXA7EmWtoIvCCpJnR9ueBCfkLycHuEkHGROBVQ865HIgz19DfonWLjyVMOvefpb54THNosGqorAy6d2/WmJxzrVODVUMKkwuNBo40sz8DnSSNyHtkJa7BqqE+fcL0Es4510Rx2gh+BYwELoq2NwB35i0iB0DXruFLf9oSwfLlXi3knMuZOIngGDP7JrAVdvUaap/XqBxSqB5KWyLwwWTOuRyKkwi2S2pL1FNIUg9gZ16jckBIBPWWCDwROOdyJE4iuAP4I7CfpCrgBeCneY3KAaGdYI8SwZYtsH69Vw0553Im06RzA8xskZlVS5oNnEroNfR/zOzdZouwhJWXw5w5KTt9MJlzLscydR99FDhK0jNmdirwXjPF5CJpSwQ+mMw5l2OZEkEbSTcAA1OXqwRfqrI5lJeHueV27IB2iX8pn2fIOZdjmdoILiT0FEqdhjrx4/KsvDxMOLduXdJOLxE453Is0zTU84H/kfSWmf21GWNykeRBZfslpvlbvhzatw9LVTrnXA7EmWvoWUkXAxXJ55vZTfkKygVpp5lYsSJUC/moYudcjsRJBH8CPgZmA5/mNxyXLO00Ez6GwDmXY3ESQV8zG533SNwe0pYIli+Hww4rSDzOudYpzoCylyQNznskbg9pp6JOVA0551yOxEkEJwCzJc2X9JaktyW9FefmkkZH130g6bp6zhklaY6kuUlrHjhCm3DXrkklgs2bw6hirxpyzuVQnKqhMY25cTQ/0Z3A6UAN8JqkqWY2L+mcfQizm442s6WSfAnMFHUGlfmoYudcHmSaYiLRP3FDI+89AvjAzBZG93sQOBeYl3TOxcDjZrYUwMxWN/K1Wq06M5D6YDLnXB5kKhHMJsw4mq6fogEHNnDv/YFlSds1wDEp5wwEyiTNIAxS+7mZ3Zt6I0kTiJbH7NevX+rhVq1HD/jww2jDB5M55/Ig04CyAU28d30JJPX1jyJMaLcX8E9JL5vZv1JimQxMBhg+fHjqPVq18nJ4881owxOBcy4P4rQRNFYNcEDSdl9geZpz1pjZJmCTpOeAocC/cMDuNQnMQCtWQIcOsO++hQ7LOdeKxOk11FivAQdLGiCpPWHuoqkp5/wJOFFSO0mdCFVHPsV1kh49YOvW0GFo1xKVPqrYOZdDeSsRmNkOSd8CngLaAlPMbK6kK6Pjk8zsXUl/A94irHp2l5m9k6+YWqLkQWV7+6hiV8S2b99OTU0NW7duLXQoJa1jx4707duXsrKy2NfE6TWUlpmty3Q8OmcaMC1l36SU7VuBWxu6V6lKnmaiYsUKH1XsilZNTQ1dunShoqICeam1IMyMtWvXUlNTw4AB8Zt5M1UNzQZmRY+1hHr796Pns5sQq8tCnWkmvETgitjWrVvp3r27J4ECkkT37t2zLpXVmwjMbICZHUio2jnHzMrNrDtwNvB4k6J1sSVKBB99uBk+/tgTgStqngQKrzH/BnEai4+OqngAiNYmOCnrV3KNkigRbF3kg8mci6OqqorDDz+cIUOGUFlZySuvvFKQOObMmcO0abtrxqdOncrNN98MwGWXXcajjz66xzUzZszg7LPPbrYYE+I0Fq+R9H+B+wnjAL4ErM1rVG6Xz30OvqRqzv/5d8OOH/wgrFs5fnxhA3OuCXr1glWr9tzfsyesXNn4+/7zn//kySef5PXXX6dDhw6sWbOGbdu2Nf6GTTBnzhxmzZrFWWedBcDYsWMZO3ZsQWJpSJwSwUVAD+CPwBPAftE+1wz0h2p+YxPosin6X7N6NUyYANXVhQ3MuSZIlwQy7Y9rxYoVlJeX06FDBwDKy8vp06cPFRUVrInmapk1axajRo0CYObMmVRWVlJZWcmwYcPYsCHMqHPLLbcwePBghg4dynXXhfkyFyxYwOjRoznqqKM48cQTee+994Dw7f7KK6/kxBNPZODAgTz55JNs27aNH/3oRzz00ENUVlby0EMPcc899/Ctb31rV6xPP/10nWtSbdq0iSuuuIKjjz6aYcOG8ac//QmAuXPnMmLECCorKxkyZAjvv/9+0940YpQIot5BVzf5lVzjTJxIJzbX3bd5M0yc6KUCV7SuuQbmzGnctdFn9B4qK+H22zNfe8YZZ3DTTTcxcOBATjvtNC644AJOOqn+muzbbruNO++8k+OPP56NGzfSsWNH/vrXv/LEE0/wyiuv0KlTJ9ZFi4ZPmDCBSZMmcfDBB/PKK6/wjW98g2effRaAxYsXM3PmTBYsWMDJJ5/MBx98wE033cSsWbP45S9/CcA999xT57XTXZOsqqqKU045hSlTprB+/XpGjBjBaaedxqRJk7j66qsZP34827Zt47PPPsv8psTQYCKQNBD4HnsuVXlKk1/dNWzp0uz2O1fCOnfuzOzZs3n++ef5xz/+wQUXXLCrXj6d448/nu985zuMHz+e8847j759+/L0009z+eWX06lTJwC6devGxo0beemllzj//PN3Xfvpp7sXbPziF79ImzZtOPjggznwwAN3lRYyaeia6dOnM3XqVG677TYg9MpaunQpI0eOpKqqipqaGs477zwOPvjgrN6jdOK0ETwCTALuApqeelx2+vWDJUvS73euSDX0zT1Tx5YZM5r22m3btmXUqFGMGjWKwYMH8/vf/5527dqxc+dOgDpdK6+77jq+8IUvMG3aNI499liefvppzGyPnjc7d+5kn332YU49xZzU8+P03GnoGjPjscceY9CgQXX2H3rooRxzzDH85S9/4cwzz+Suu+7ilFOa9r08ThvBDjP7tZm9amazEz9NelUXX1UVn7brVHdfp05QVVWYeJwrYvPnz69TZz5nzhz69+9PRUUFs2eHj63HHnts1/EFCxYwePBgrr32WoYPH857773HGWecwZQpU9i8OVTJrlu3jq5duzJgwAAeeeQRIHxIv7lrNkh45JFH2LlzJwsWLGDhwoUMGjSILl267GpzSCfdNcnOPPNMfvGLX2AW5tl84403AFi4cCEHHnggV111FWPHjuWtt2KtE5ZRnETwZ0nfkNRbUrfET5Nf2cUzfjx/PmcyW+kQpm7t3x8mT/b2Adei9eyZ3f64Nm7cyKWXXsphhx3GkCFDmDdvHjfeeCM33HADV199NSeeeCJt27bddf7tt9/OEUccwdChQ9lrr70YM2YMo0ePZuzYsQwfPpzKyspdVTPV1dX87ne/Y+jQoRx++OG7Gm8BBg0axEknncSYMWOYNGkSHTt25OSTT2bevHm7GotTpbsm2fXXX8/27dsZMmQIRxxxBNdffz0ADz30EEcccQSVlZW89957XHLJJU170wAlsk29J0iL0uy2aLBZsxs+fLjNmjWrEC9dML/8JZzx7YH0+/ej6Pj4A4UOx7m03n33XQ499NBCh9HsLrvsMs4++2zGjRtX6FB2SfdvIWm2mQ1Pd36cXkNNXZfANVF5OZSzhs17ldOx4dOdcy4rsWYflXQEcBjs/hxKt5KYy4/99t1ONz5icYceeJ2cc8UltVtoSxSn++gNwChCIphGWMz+BcATQTPpWRb6Ma9vV17gSJxzrVGcxuJxhKUkV5rZ5YQVxDrkNSpXRw9qAag1TwTOudyLkwi2mNlOYIekrsBqGl643uXQPjvC0PhVO3sUOBLnXGsUp41glqR9gN8S1iHYCLyaz6BcXe0/DiWCDz/1EoFzLvcaLBGY2TfMbH20stjpwKVRFZFrLtFkWUu3eInAuUw6d+7c4DnPP/88hx9+OJWVlWzZsiWr+z/xxBPMmzdv1/aPfvQjnn766azjLDZZLV5vZovNrOnD2Fx2okSweEP3AgfiXA5VV0NFBbRpEx6baUbd6upqvve97zFnzhz22muvrK5NTQQ33XQTp512Wq5DbHZZJQLX/Hr1gjt+VMtH7MO0v5chhXlaevUqdGTONUF1dZhOfckSMAuPOZxefcaMGYwaNYpx48ZxyCGHMH78eMyMu+66i4cffpibbrqJ8dHo/FtvvZWjjz6aIUOGcMMNN+y6x7333suQIUMYOnQoX/7yl3nppZeYOnUq3//+96msrGTBggV1Fph55plnGDZsGIMHD+aKK67YNSldtlNgF0KscQSucFatCoPJ1lC+x37nilZD81C//DIkzd4JhOnVv/IV+O1v018TZx7qJG+88QZz586lT58+HH/88bz44ot89atf5YUXXtg1Enj69Om8//77vPrqq5gZY8eO5bnnnqN79+5UVVXx4osvUl5ezrp16+jWrRtjx45NO4p469atXHbZZTzzzDMMHDiQSy65hF//+tdcc8019caXbgrsQmmwRCDpIEkdouejJF0VNR67ZtKDWmrx9gHXiqQmgYb2N8KIESPo27cvbdq0obKyksWLF+9xzvTp05k+fTrDhg3jyCOP5L333uP999/n2WefZdy4cZRHa8V265Z5KOf8+fMZMGAAAwcOBODSSy/lueeey3hNYgrsO+64g/Xr19OuXeG+l8d55ceA4ZL+DfgdMBX4A3BWPgNzu5WzhmUcUOgwnIuvoW/uFRXpp1fv37/p81BHEquUQZiaeseOHXucY2b88Ic/5Gtf+1qd/XfccUdWi8BnmrMtmymwDznkkNivmUtx2gh2mtkO4N+B283sPwFfQb0Z9aB2j6oh51q0qqownXqyAkyvfuaZZzJlyhQ2btwIwIcffsjq1as59dRTefjhh1m7NizPnlilrL6ppQ855BAWL168a5Wx++67b9fKaNlMgV0ocRLBdkkXAZcCiYU1y/IXkqvLKGeNVw251mX8+DCdev/+ofdDgaZXP+OMM7j44osZOXIkgwcPZty4cWzYsIHDDz+ciRMnctJJJzF06FC+853vAHDhhRdy6623MmzYMBYsWLDrPh07duTuu+/m/PPPZ/DgwbRp04Yrr7wSIKspsAslzjTUhwFXAv80swckDQAuMLP613/Lo1Kbhvqg/TawoLYr3+cWbuP7u/b37AkrVxYwMOdSlOo01MUo22mo4wwomwdcC7webS8qVBIoRQteCd3ObvpVKBHcemvobedJwDmXK3F6DZ0DzAH+Fm1XSpqa57hcQm2YXqJj33L23huWLy9wPM65VidOG8GNwAhgPYCZzQF8sZrmEg1E0X496NMHVqwocDzOuVYn7uL1H6fsy9yw4HInSgSUl9O7t5cIXHFrqM3R5V9j/g3iJIJ3JF0MtJV0sKRfAC9l/UqucaKqIXqEEoEnAlesOnbsyNq1az0ZFJCZsXbt2qxHKccZUPZtYCLwKfAA8BTwkzg3lzQa+DnQFrirvkZmSUcDLxN6Iz0a594lY80aKCuDLl12VQ2ZhR53zhWTvn37UlNTQ23iy4sriI4dO9K3b9+sromzeP1mQiKYmM2NJbUF7iRMXV0DvCZpatQLKfW8/yEkGJeqtjasXi/Ruzds2gQbNkDXroUOzLm6ysrKGDDAmw9bonoTgaQ/k6EtwMzGNnDvEcAHZrYwut+DwLnAvJTzvk2YxuLoOAGXnDVroEfoOtqnT9i1fLknAudc7mQqEdzWxHvvDyxL2q4Bjkk+QdL+hKkrTiFDIpA0AZgA0K9fvyaG1cIkSgTsTgQrVkCBpiRxzrVC9SYCM5uZeC6pPXAIoYQw38y2xbh3ulrs1BLG7cC1ZvZZpgmezGwyMBnCyOIYr916rFkDw4YB0Dua4ckbjJ1zudRgG4GkLwCTgAWED/cBkr5mZn9t4NIaqDNlZl8g9SNsOPBglATKgbMk7TCzJ+KFXwLSlAg8ETjncilOr6GfASeb2QcQ1icA/gI0lAheAw6O5ib6ELgQuDj5BDPb1bIk6R7gSU8CSXbsgI8+2tVG0KUL7L23DypzzuVWnESwOpEEIguB1Q1dZGY7JH2L0BuoLTDFzOZKujI6PqkxAZeUaOrbRIkA8LEEzrmcy9Rr6Lzo6VxJ04CHCXX85xO+7TfIzKYB01L2pU0AZnZZnHuWlKTBZAmeCJxzuZapRHBO0vNVwEnR81pg37xF5HZLml4ioXdvKKFZuJ1zzSBTr6HLmzMQl0aiRJCmashHFzvnciVT1dAPzOyWaG6hPbpsmtlVeY3M7S4RpFQNbd4Mn3wCn/tcgeJyzrUqmaqG3o0evSKiUBIlgu7dd+1KHkvgicA5lwuZqob+HD3dbGaPJB+TdH5eo3LBmjXh0759+127kkcX+6qAzrlciDMN9Q9j7nO5ljSYLMEHlTnnci1TG8EY4Cxgf0l3JB3qCuzId2COOhPOJfg0E865XMvURrCc0D4wFpidtH8D8J/5DMpFamshZV7xLl2gc2cfXeycy51MbQRvAm9K+gNhjqFsJ51zTbVmDVRW7rHbB5U553IpzhQTpwO/IftJ51xTmKWtGgJ87WLnXE7FSQT/S+MmnXNNsWkTbN26R2MxhBLBq68WICbnXKsUp9dQoyadc02UZjBZQvLoYueca6o4JYK0k84lJqUzs8fzGF/pSjO9RELv3rBlC3z8MeyzT/OG5ZxrfeIkgo7sOelcN8KkdAZ4IsiHBkoEEHoOeSJwzjVVg4nAJ58rkAwlguRBZT662DnXVA22EUgaKOkZSe9E20Mk/d/8h1biMpQIfFCZcy6X4jQW/5YwpcR2ADN7i7DspMun2lpo1w66dt3jUCIR+KAy51wuxEkEncwstbOiTzGRb2vWhGqhNIsOdOkSfrxE4JzLhTiJYE00dsAAJI0D/LtoPlVXh5+VK6GiIjxP4YPKnHO5EqfX0DeBycAhkj4EFgFfymtUpay6GiZMCP1DAZYsCdsA48fvOq1PH68acs7lRoMlAjNbaGanAT2AQ8zsBDNbnPfIStXEiWEJsmSbN4f9SXy+IedcrsTpNfRTSfuY2SYz2yBpX0n/3RzBlaSlS2PtT1QN+ehi51xTxWkjGGNm6xMbZvYRYZ0Clw/9+jW4v1cv+NnPwlREbdqE9mQp7HfOuWzFSQRtJXVIbEjaC+iQ4XzXFFVV0CHl7e3UKeyPrFqV/tL69jvnXCZxEsH9wDOSviLpCuDvwO/zG1YJGz8errgiPJegf3+YPLlOQ7FzzuVSnCkmbpH0FnAaYT2Cn5jZU3mPrJQl5pDYtAn22quwsTjnWr0GE4GkvYHpZvY3SYOAQZLKzGx7/sMrUYsWhQp/TwLOuWYQp2roOaCjpP2Bp4HLgXvyGVTJW7QIBgwodBTOuRIRJxHIzDYD5wG/MLN/Bw7Lb1glroFE0LNndvudcy6TWIlA0khgPGGJSog3Itk1xo4dsGxZxkSwcmUYP/Dii2H7j38M2ytXNlOMzrlWJU4iuIYw++gfzWyupAOBf8S5uaTRkuZL+kDSdWmOj5f0VvTzkqShWUXfGi1bBp99Fqtq6MgjoX373QnBOecaI06voZnAzKTthcBVDV0nqS1wJ3A6UENY3nKqmc1LOm0RcJKZfSRpDGFOo2Oy+xVamcWLw2OMRNCxIwwfDi+9lN+QnHOtW5xeQ/8gmnk0mZmd0sClI4APosSBpAeBc4FdicDMkj/CXgb6xoi5dVu0KDxWVMQ6/bjj4I474NNP9xyH5pxzccSpGvoe8P3o53pgDjArxnX7A8uStmuiffX5CvDXdAckTZA0S9Ks2sQSjq3VokVh3ogDDoh1+nHHwbZtMHt2nuNyzrVacaqGUj9iXpQ0M+3Jde25okqakgWApJMJieCEemKYTKg2Yvjw4a17mrVFi0ISKCuLdfpxx4XHl17a/dw557IRp2qoW9JmG+AoIM70ZjVA8tfavsAeEydLGgLcRZjcbm2M+7ZuWY4h6NkTDjrI2wmcc40XpxvobMI3eRGWqFxE+PbekNeAgyUNAD4krHN8cfIJkvoBjwNfNrN/ZRF367VoEYwendUlxx0HTz0VupCmWdnSOecyilM11Kghrma2Q9K3gKeAtsCUqPvpldHxScCPgO7ArxQ+wXaY2fDGvF6rsGVLWHYsy1HFxx8P990HCxeG0oFzzmWj3kQg6WhgmZmtjLYvAf4DWALcaGbrGrq5mU0DpqXsm5T0/KvAVxsXeiu0ZEl4zDIRJLcTeCJwzmUrU6+h3wDbACR9HrgZuBf4mKjh1uVYFmMIkh12GHTt6gPLnHONk6lqqG3St/4LgMlm9hjwmKQ5eY+sFCXGEGSZCNq2hZEjvcHYOdc4GROBpHZmtgM4FZgQ8zrXWIsWhVFhWa452avX7tXJkhuLe/b0+Yeccw3L9IH+ADBT0hpgC/A8gKR/I1QPuVxbtCisSNYmzji/3XzpSudcU9SbCMysStIzQG/CwjSJgVxtgG83R3Alx9chcM4VQMYqHjN7Oc0+7++fL4sWwdFHFzoK51yJya4OwuXPJ5/AunVeInDONTtPBMWikT2GnHOuqTwRFItGjiGA+peo7NGj8eE450qHJ4Ji0YQSQWLpysRPeXnYX1sbupMmfrLsleqcKxGeCIrFokXQpQt069bwuQ1Ysyb9/lWr6iYGTw7OOfCBYcVj0aKwKlkzTx+aSA7OuZYllwNGvURQDKqr4W9/g7ffDsmgurrQETnnilwuB4x6Iii06mqYMAG2bw/bS5aEbU8Gzrlm4omg0CZOhM2b6+7bvDnsd865ZuCJoNCWLs1ufwz1dSd1zrl0PBEUWr9+2e2PIbU7qZknB+dc/TwRFNpPfrLnvk6doKoqpy+Tmhw8MTjXsuXy/7B3Hy20vn3DY3k5rF0bSgJVVTB+fF5f1tcpcM4leCIotIceCiWAJUvCo3PONTOvGiqkHTvgscfgnHM8CTjnCsYTQSHNmBHmg7jggkJH4pwrYZ4ICumhh6BzZxg9utCROOdKmCeCQtm+HR5/HMaOhb32KnQ0zrkS5omgUJ59NqxI5tVCzrkC80TQ3Kqrw8Ryo0eHaT/XrSt0RM65EufdR5tTYoK5xNxCZvDNb0JZWd7HDTjnXH28RNCcfII551wR8kTQnPIwwZxzzjWVJ4J8S7QJtGlT/1JgTZhgzjnnmsoTQa4lf/CXl8MVV4TpI8xg5849z8/DBHPOOZeNvCYCSaMlzZf0gaTr0hyXpDui429JOjIvgSR/OCeWgkzd941vZN6Oe82ECbs/+NeuhW3b9oynbdtQOujfHyZP9oZi51xByczyc2OpLfAv4HSgBngNuMjM5iWdcxbwbeAs4Bjg52Z2TKb7Dh8+3GbNmhU/kNSeOhB66UjpP6TrE+caKSSAhkjpSwfOOZcnkmab2fB0x/JZIhgBfGBmC81sG/AgcG7KOecC91rwMrCPpN45jSJdT53t27NLAnGviZtUvU3AOVdE8pkI9geWJW3XRPuyPQdJEyTNkjSrtrY2uyiKrUeOtwk454pMPhNBui4yqV+Z45yDmU02s+FmNrxHjx7ZRdHc375TewaVlUH37t4m4JwrWvlMBDXAAUnbfYHljTinaaqq9pzrv6wM2rfP7j5xrunUCa68MnzgJz747747TDW9cycsXuxJwDlXdPKZCF4DDpY0QFJ74EJgaso5U4FLot5DxwIfm9mKnEYxfnz4Fp764TxlSt19X/965u0410yeDL/6VfjA9w9+51wLkbdeQ7CrV9DtQFtgiplVSboSwMwmSRLwS2A0sBm43MwydgnKuteQc865jL2G8jrpnJlNA6al7JuU9NyAb+YzBuecc5n5yGLnnCtxngicc67EeSJwzrkS54nAOedKXF57DeWDpFpgSSMvLwfW5DCcfGtJ8bakWKFlxduSYoWWFW9LihWaFm9/M0s7IrfFJYKmkDSrvu5TxaglxduSYoWWFW9LihVaVrwtKVbIX7xeNeSccyXOE4FzzpW4UksEkwsdQJZaUrwtKVZoWfG2pFihZcXbkmKFPMVbUm0Ezjnn9lRqJQLnnHMpPBE451yJK5lEIGm0pPmSPpB0XaHjSSVpiqTVkt5J2tdN0t8lvR897lvIGBMkHSDpH5LelTRX0tXR/qKLV1JHSa9KejOK9cfFGmuCpLaS3pD0ZLRdzLEulvS2pDmSZkX7ijnefSQ9Kum96O93ZDHGK2lQ9J4mfj6RdE2+Yi2JRCCpLXAnMAY4DLhI0mGFjWoP9xCm4052HfCMmR0MPBNtF4MdwHfN7FDgWOCb0ftZjPF+CpxiZkOBSmB0tPZFMcaacDXwbtJ2MccKcLKZVSb1by/meH8O/M3MDgGGEt7noovXzOZH72klcBRhmv4/kq9YzazV/wAjgaeStn8I/LDQcaWJswJ4J2l7PtA7et4bmF/oGOuJ+0/A6cUeL9AJeB04plhjJazS9wxwCvBksf8dAIuB8pR9RRkv0BVYRNRJptjjTYrvDODFfMZaEiUCYH9gWdJ2TbSv2PW0aMW26HG/AsezB0kVwDDgFYo03qiqZQ6wGvi7mRVtrISFnH4A7EzaV6yxQlhjfLqk2ZImRPuKNd4DgVrg7qjq7S5Je1O88SZcCDwQPc9LrKWSCJRmn/ebbSJJnYHHgGvM7JNCx1MfM/vMQhG7LzBC0hEFDiktSWcDq81sdqFjycLxZnYkodr1m5I+X+iAMmgHHAn82syGAZsogmqgTKJlfscCj+TzdUolEdQAByRt9wWWFyiWbKyS1Bsgelxd4Hh2kVRGSALVZvZ4tLto4wUws/XADEJbTDHGejwwVtJi4EHgFEn3U5yxAmBmy6PH1YQ67BEUb7w1QE1UIgR4lJAYijVeCAn2dTNbFW3nJdZSSQSvAQdLGhBl2AuBqQWOKY6pwKXR80sJdfEFF601/TvgXTP736RDRRevpB6S9ome7wWcBrxHEcZqZj80s75mVkH4G33WzL5EEcYKIGlvSV0Szwl12e9QpPGa2UpgmaRB0a5TgXkUabyRi9hdLQT5irXQDSHN2OByFvAvYAEwsdDxpInvAWAFsJ3wzeUrQHdCw+H70WO3QscZxXoCoWrtLWBO9HNWMcYLDAHeiGJ9B/hRtL/oYk2JexS7G4uLMlZCnfub0c/cxP+rYo03iq0SmBX9PTwB7Fus8RI6N6wFPpe0Ly+x+hQTzjlX4kqlasg551w9PBE451yJ80TgnHMlzhOBc86VOE8EzjlX4jwROFcPSd2TZn9cKenD6PlGSb8qdHzO5Yp3H3UuBkk3AhvN7LZCx+JcrnmJwLksSRqVtFbAjZJ+L2l6NDf/eZJuiebo/1s0FQeSjpI0M5qc7anENAHOFQNPBM413UHAF4BzgfuBf5jZYGAL8IUoGfwCGGdmRwFTgKpCBetcqnaFDsC5VuCvZrZd0ttAW+Bv0f63CWtMDAKOAP4epmmiLWE6EeeKgicC55ruUwAz2ylpu+1ueNtJ+D8mYK6ZjSxUgM5l4lVDzuXffKCHpJEQpvCWdHiBY3JuF08EzuWZmW0DxgH/I+lNwmytxxU0KOeSePdR55wrcV4icM65EueJwDnnSpwnAuecK3GeCJxzrsR5InDOuRLnicA550qcJwLnnCtx/x/Q1DY3zF2Y6QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "\n",
    "# p2.1-SIS\n",
    "# Second, the python script for SIS\n",
    "\n",
    "# -*- coding: utf-8 -*-\n",
    "\n",
    "import scipy.integrate as spi\n",
    "import numpy as np\n",
    "import pylab as pl\n",
    "\n",
    "beta=1.4247\n",
    "gamma=0.14286\n",
    "I0=1e-6\n",
    "ND=70\n",
    "TS=1.0\n",
    "INPUT = (1.0-I0, I0)\n",
    "\n",
    "def diff_eqs(INP,t):\n",
    "\t'''The main set of equations'''\n",
    "\tY=np.zeros((2))\n",
    "\tV = INP\n",
    "\tY[0] = - beta * V[0] * V[1] + gamma * V[1]\n",
    "\tY[1] = beta * V[0] * V[1] - gamma * V[1]\n",
    "\treturn Y   # For odeint\n",
    "\n",
    "t_start = 0.0; t_end = ND; t_inc = TS\n",
    "t_range = np.arange(t_start, t_end+t_inc, t_inc)\n",
    "RES2 = spi.odeint(diff_eqs,INPUT,t_range)\n",
    "\n",
    "\n",
    "#Ploting\n",
    "pl.plot(RES2[:,0], '-bs', label='Susceptibles')\n",
    "pl.plot(RES2[:,1], '-ro', label='Infectious')\n",
    "pl.legend(loc=0)\n",
    "pl.title('SIS epidemic without births or deaths')\n",
    "pl.xlabel('Time')\n",
    "pl.ylabel('Susceptibles and Infectious')\n",
    "#pl.savefig('2.5-SIS-high.png', dpi=900) # This does increase the resolution.\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-23T12:57:27.330597Z",
     "start_time": "2021-09-23T12:57:27.311020Z"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df = pd.DataFrame(RES2, columns = ['Susceptibles', 'Infectious'])\n",
    "df.to_excel('Fig2.5-SIS.xlsx', index = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
